← Back to branch summary

~ubuntu-branches/ubuntu/oneiric/mpqc/oneiric

« back to all changes in this revision

Viewing changes to src/lib/chemistry/qc/wfn/wfn.cc

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck
  • Date: 2005-11-27 11:41:49 UTC
  • mfrom: (1.1.2 upstream)
  • Revision ID: james.westby@ubuntu.com-20051127114149-zgz9r3gk50w8ww2q
Tags: 2.3.0-1
http://bugs.debian.org/339232
* New upstream release.
* debian/rules (SONAME): Activate awk snippet for automatic so-name
  detection again, resulting in a bump to `7' and making a `c2a' for
  the C++ allocator change unnecessary; closes: #339232.
* debian/patches/00list (08_gcc-4.0_fixes): Removed, no longer needed.
* debian/rules (test): Remove workarounds, do not abort build if tests
  fail.
* debian/ref: Removed.
* debian/control.in (libsc): Added Conflict against libsc6c2.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/lib/chemistry/qc/wfn/wfn.cc
32
32
#include <stdlib.h>
33
33
#include <math.h>
34
34
 
 
35
#include <stdexcept>
35
36
#include <iostream>
 
37
#include <iterator>
36
38
 
37
39
#include <util/keyval/keyval.h>
38
40
#include <util/misc/timer.h>
39
41
#include <util/misc/formio.h>
 
42
#include <util/misc/autovec.h>
40
43
#include <util/state/stateio.h>
41
44
#include <chemistry/qc/basis/obint.h>
42
45
#include <chemistry/qc/basis/symmint.h>
49
52
 
50
53
#define CHECK_SYMMETRIZED_INTEGRALS 0
51
54
 
 
55
/////////////////////////////////////////////////////////////////////////
 
56
 
 
57
// This maps a TwoBodyThreeCenterInt into a OneBodyInt
 
58
namespace sc {
 
59
class ChargeDistInt: public OneBodyInt {
 
60
    Ref<TwoBodyThreeCenterInt> tbint_;
 
61
    Ref<Molecule> mol_;
 
62
    Ref<GaussianBasisSet> ebasis0_;
 
63
    Ref<GaussianBasisSet> ebasis1_;
 
64
    Ref<GaussianBasisSet> atom_basis_;
 
65
    std::vector<int> i_cd_;
 
66
    const double *coef_;
 
67
  public:
 
68
    // The electronic basis are bs1 and bs2 in tbint and the
 
69
    // nuclear basis is bs3.
 
70
    ChargeDistInt(const Ref<TwoBodyThreeCenterInt>& tbint,
 
71
                  const double *coef);
 
72
 
 
73
    void compute_shell(int,int);
 
74
 
 
75
    bool cloneable();
 
76
};
 
77
 
 
78
ChargeDistInt::ChargeDistInt(const Ref<TwoBodyThreeCenterInt>& tbint,
 
79
                             const double *coef):
 
80
  OneBodyInt(tbint->integral(),tbint->basis1(),tbint->basis2()),
 
81
  tbint_(tbint),
 
82
  coef_(coef)
 
83
{
 
84
  ebasis0_ = tbint->basis1();
 
85
  ebasis1_ = tbint->basis2();
 
86
  atom_basis_ = tbint->basis3();
 
87
  mol_ = atom_basis_->molecule();
 
88
  buffer_ = new double[tbint->basis1()->max_nfunction_in_shell()
 
89
                       *tbint->basis2()->max_nfunction_in_shell()];
 
90
 
 
91
  for (int i=0; i<atom_basis_->ncenter(); i++) {
 
92
    if (atom_basis_->nshell_on_center(i) > 0) i_cd_.push_back(i);
 
93
  }
 
94
}
 
95
 
 
96
void
 
97
ChargeDistInt::compute_shell(int ish,int jsh)
 
98
{
 
99
//   std::cout << "starting " << ish << " " << jsh << std::endl;
 
100
  int nijbf
 
101
    = ebasis0_->shell(ish).nfunction()
 
102
    * ebasis1_->shell(jsh).nfunction();
 
103
  memset(buffer_,0,nijbf*sizeof(buffer_[0]));
 
104
  const double *tbbuffer = tbint_->buffer();
 
105
  int ksh = 0;
 
106
  int coef_offset = 0;
 
107
  for (int ii=0; ii<i_cd_.size(); ii++) {
 
108
    int i = i_cd_[ii];
 
109
    int nshell = atom_basis_->nshell_on_center(i);
 
110
    for (int j=0; j<nshell; j++, ksh++) {
 
111
      std::cout << "computing " << ish << " " << jsh << " " << ksh << std::endl;
 
112
      tbint_->compute_shell(ish,jsh,ksh);
 
113
      int nbfk = atom_basis_->shell(i).nfunction();
 
114
      for (int ijbf=0; ijbf<nijbf; ijbf++) {
 
115
        for (int kbf=0; kbf<nbfk; kbf++) {
 
116
          buffer_[ijbf]
 
117
            -= coef_[coef_offset+kbf] * tbbuffer[ijbf*nbfk + kbf];
 
118
//           std::cout << "  adding "
 
119
//                     << coef_[coef_offset+kbf] * tbbuffer[ijbf*nbfk + kbf]
 
120
//                     << " = "
 
121
//                     << coef_[coef_offset+kbf]
 
122
//                     << " * "
 
123
//                     << tbbuffer[ijbf*nbfk + kbf]
 
124
//                     << " at location "
 
125
//                     << ijbf
 
126
//                     << std::endl;
 
127
        }
 
128
      }
 
129
    }
 
130
    coef_offset += nshell;
 
131
  }
 
132
//   std::cout << "done with " << ish << " " << jsh << std::endl;
 
133
}
 
134
 
 
135
bool
 
136
ChargeDistInt::cloneable()
 
137
{
 
138
  // not cloneable because tbint is not cloneable
 
139
  return false;
 
140
}
 
141
 
 
142
}
 
143
 
 
144
/////////////////////////////////////////////////////////////////////////
 
145
 
52
146
static ClassDesc Wavefunction_cd(
53
 
  typeid(Wavefunction),"Wavefunction",5,"public MolecularEnergy",
 
147
  typeid(Wavefunction),"Wavefunction",7,"public MolecularEnergy",
54
148
  0, 0, 0);
55
149
 
56
150
Wavefunction::Wavefunction(const Ref<KeyVal>&keyval):
119
213
    "Wavefunction::Wavefunction\n"
120
214
    );
121
215
 
 
216
  atom_basis_ << keyval->describedclassvalue("atom_basis");
 
217
  if (atom_basis_.nonnull()) {
 
218
    atom_basis_coef_ = new double[atom_basis_->nbasis()];
 
219
    for (int i=0; i<atom_basis_->nbasis(); i++) {
 
220
      atom_basis_coef_[i] = keyval->doublevalue("atom_basis_coef",i);
 
221
    }
 
222
    scale_atom_basis_coef();
 
223
 
 
224
  }
 
225
  else {
 
226
    atom_basis_coef_ = 0;
 
227
  }
 
228
 
122
229
  integral_ << keyval->describedclassvalue("integrals");
123
230
  if (integral_.null()) {
124
231
    Integral* default_intf = Integral::get_default_integral();
189
296
  gbs_ << SavableState::restore_state(s);
190
297
  integral_ << SavableState::restore_state(s);
191
298
 
 
299
  if (s.version(::class_desc<Wavefunction>()) >= 6) {
 
300
    Ref<KeyVal> original_override = s.override();
 
301
    Ref<AssignedKeyVal> matrixkit_override = new AssignedKeyVal;
 
302
    matrixkit_override->assign("matrixkit", gbs_->so_matrixkit().pointer());
 
303
    Ref<KeyVal> new_override
 
304
      = new AggregateKeyVal(matrixkit_override.pointer(),
 
305
                            original_override.pointer());
 
306
    s.set_override(new_override);
 
307
    orthog_ << SavableState::restore_state(s);
 
308
    s.set_override(original_override);
 
309
  }
 
310
 
 
311
  if (s.version(::class_desc<Wavefunction>()) >= 7) {
 
312
    atom_basis_ << SavableState::restore_state(s);
 
313
    if (atom_basis_.nonnull()) {
 
314
      s.get_array_double(atom_basis_coef_, atom_basis_->nbasis());
 
315
    }
 
316
  }
 
317
  else {
 
318
    atom_basis_coef_ = 0;
 
319
  }
 
320
 
192
321
  integral_->set_basis(gbs_);
193
322
  Ref<PetiteList> pl = integral_->petite_list();
194
323
 
238
367
 
239
368
  SavableState::save_state(gbs_.pointer(),s);
240
369
  SavableState::save_state(integral_.pointer(),s);
 
370
  SavableState::save_state(orthog_.pointer(),s);
 
371
  SavableState::save_state(atom_basis_.pointer(), s);
 
372
  if (atom_basis_.nonnull()) {
 
373
    s.put_array_double(atom_basis_coef_,atom_basis_->nbasis());
 
374
  }
241
375
}
242
376
 
243
377
double
263
397
 
264
398
      RefSymmSCMatrix densortho(oso_dimension(), basis_matrixkit());
265
399
      densortho.assign(0.0);
266
 
      densortho.accumulate_transform(so_to_orthog_so(),dens);
 
400
      densortho.accumulate_transform(so_to_orthog_so_inverse().t(),dens);
267
401
 
268
402
      RefSCMatrix natorb(oso_dimension(), oso_dimension(),
269
403
                         basis_matrixkit());
385
519
    hao.element_op(hc);
386
520
    hc=0;
387
521
 
388
 
    Ref<OneBodyInt> nuc = integral_->nuclear();
389
 
    nuc->reinitialize();
390
 
    hc = new OneBodyIntOp(new SymmOneBodyIntIter(nuc, pl));
391
 
    hao.element_op(hc);
392
 
    hc=0;
 
522
    if (atom_basis_.null()) {
 
523
      Ref<OneBodyInt> nuc = integral_->nuclear();
 
524
      nuc->reinitialize();
 
525
      hc = new OneBodyIntOp(new SymmOneBodyIntIter(nuc, pl));
 
526
      hao.element_op(hc);
 
527
      hc=0;
 
528
    }
 
529
    else {
 
530
      // we have an atom_basis, so some nuclear charges will be computed
 
531
      // with the two electron integral code and some will be computed
 
532
      // with the point charge code
 
533
 
 
534
      // find out which atoms are point charges and which ones are charge
 
535
      // distributions
 
536
      std::vector<int> i_pc; // point charge centers
 
537
      for (int i=0; i<atom_basis_->ncenter(); i++) {
 
538
        if (atom_basis_->nshell_on_center(i) == 0) i_pc.push_back(i);
 
539
      }
 
540
      int n_pc = i_pc.size();
 
541
 
 
542
      // initialize the point charge data
 
543
      auto_vec<double> pc_charges(new double[n_pc]);
 
544
      auto_vec<double*> pc_xyz(new double*[n_pc]);
 
545
      auto_vec<double> pc_xyz0(new double[n_pc*3]);
 
546
      pc_xyz[0] = pc_xyz0.get();
 
547
      for (int i=1; i<n_pc; i++) pc_xyz[i] = pc_xyz[i-1] + 3;
 
548
      for (int i=0; i<n_pc; i++) {
 
549
        pc_charges[i] = molecule()->charge(i_pc[i]);
 
550
        for (int j=0; j<3; j++) pc_xyz[i][j] = molecule()->r(i_pc[i],j);
 
551
      }
 
552
      Ref<PointChargeData> pc_data
 
553
        = new PointChargeData(n_pc, pc_xyz.get(), pc_charges.get());
 
554
 
 
555
      // compute the point charge contributions
 
556
      Ref<OneBodyInt> pc_int = integral_->point_charge(pc_data);
 
557
      hc = new OneBodyIntOp(new SymmOneBodyIntIter(pc_int,pl));
 
558
      hao.element_op(hc);
 
559
      hc=0;
 
560
      pc_int=0;
 
561
 
 
562
      // compute the charge distribution contributions
 
563
      // H_ij += sum_A -q_A sum_k N_A_k (ij|A_k)
 
564
      integral()->set_basis(gbs_,gbs_,atom_basis_);
 
565
      Ref<TwoBodyThreeCenterInt> cd_tbint
 
566
        = integral_->electron_repulsion3();
 
567
      Ref<OneBodyInt> cd_int = new ChargeDistInt(cd_tbint, atom_basis_coef_);
 
568
      hc = new OneBodyIntOp(new SymmOneBodyIntIter(cd_int,pl));
 
569
      hao.element_op(hc);
 
570
      integral()->set_basis(gbs_);
 
571
      hc=0;
 
572
      cd_int=0;
 
573
      cd_tbint=0;
 
574
    }
393
575
 
394
576
    // now symmetrize Hso
395
577
    RefSymmSCMatrix h(so_dimension(), basis_matrixkit());
407
589
    hao.element_op(hc);
408
590
    hc=0;
409
591
 
 
592
//     std::cout << "null atom_basis" << std::endl;
410
593
    Ref<OneBodyInt> nuc = integral_->nuclear();
411
594
    nuc->reinitialize();
412
595
    hc = new OneBodyIntOp(new OneBodyIntIter(nuc));
459
642
  return hcore_.result_noupdate();
460
643
}
461
644
 
 
645
// Compute lists of centers that are point charges and lists that
 
646
// are charge distributions.
 
647
void
 
648
Wavefunction::set_up_charge_types(
 
649
  std::vector<int> &q_pc,
 
650
  std::vector<int> &q_cd,
 
651
  std::vector<int> &n_pc,
 
652
  std::vector<int> &n_cd)
 
653
{
 
654
  q_pc.clear();
 
655
  q_cd.clear();
 
656
  n_pc.clear();
 
657
  n_cd.clear();
 
658
 
 
659
  for (int i=0; i<atom_basis_->ncenter(); i++) {
 
660
    bool is_Q = (molecule()->atom_symbol(i) == "Q");
 
661
    if (atom_basis_->nshell_on_center(i) > 0) {
 
662
      if (is_Q) q_cd.push_back(i);
 
663
      else      n_cd.push_back(i);
 
664
    }
 
665
    else {
 
666
      if (is_Q) q_pc.push_back(i);
 
667
      else      n_pc.push_back(i);
 
668
    }
 
669
  }
 
670
}
 
671
 
 
672
double
 
673
Wavefunction::nuclear_repulsion_energy()
 
674
{
 
675
  if (atom_basis_.null()) return molecule()->nuclear_repulsion_energy();
 
676
 
 
677
  double nucrep = 0.0;
 
678
 
 
679
  std::vector<int> q_pc, q_cd, n_pc, n_cd;
 
680
  set_up_charge_types(q_pc,q_cd,n_pc,n_cd);
 
681
 
 
682
  // compute point charge - point charge contribution
 
683
  nucrep += nuc_rep_pc_pc(n_pc, n_pc, true /* i > j  */);
 
684
  nucrep += nuc_rep_pc_pc(q_pc, n_pc, false /* all i j */);
 
685
  if (molecule()->include_qq()) {
 
686
    nucrep += nuc_rep_pc_pc(q_pc, q_pc, true /* i > j */);
 
687
  }
 
688
 
 
689
  // compute point charge - charge distribution contribution
 
690
  nucrep += nuc_rep_pc_cd(n_pc, n_cd);
 
691
  nucrep += nuc_rep_pc_cd(q_pc, n_cd);
 
692
  nucrep += nuc_rep_pc_cd(n_pc, q_cd);
 
693
  if (molecule()->include_qq()) {
 
694
    nucrep += nuc_rep_pc_cd(q_pc, q_cd);
 
695
  }
 
696
 
 
697
  // compute the charge distribution - charge distribution contribution
 
698
  nucrep += nuc_rep_cd_cd(n_cd, n_cd, true /* i > j  */);
 
699
  nucrep += nuc_rep_cd_cd(q_cd, n_cd, false /* all i j  */);
 
700
  if (molecule()->include_qq()) {
 
701
    nucrep += nuc_rep_cd_cd(q_cd, q_cd, true /* i > j */);
 
702
  }
 
703
 
 
704
  return nucrep;
 
705
}
 
706
 
 
707
double
 
708
Wavefunction::nuc_rep_pc_pc(const std::vector<int>&c1,
 
709
                            const std::vector<int>&c2,
 
710
                            bool uniq)
 
711
{
 
712
  double e = 0.0;
 
713
 
 
714
  if (c1.size() == 0 || c2.size() == 0) return e;
 
715
 
 
716
  for (int ii=0; ii<c1.size(); ii++) {
 
717
    int i = c1[ii];
 
718
    SCVector3 ai(molecule()->r(i));
 
719
    double Zi = molecule()->charge(i);
 
720
    int jfence = (uniq?ii:c2.size());
 
721
    for (int jj=0; jj<jfence; jj++) {
 
722
      int j = c2[jj];
 
723
      SCVector3 aj(molecule()->r(j));
 
724
      e += Zi * molecule()->charge(j) / ai.dist(aj);
 
725
    }
 
726
  }
 
727
 
 
728
  return e;
 
729
}
 
730
 
 
731
double
 
732
Wavefunction::nuc_rep_pc_cd(const std::vector<int>&pc,
 
733
                            const std::vector<int>&cd)
 
734
{
 
735
  double e = 0.0;
 
736
 
 
737
  if (pc.size() == 0 || cd.size() == 0) return e;
 
738
 
 
739
  integral()->set_basis(atom_basis());
 
740
 
 
741
  sc::auto_vec<double> charges(new double[pc.size()]);
 
742
  sc::auto_vec<double*> positions(new double*[pc.size()]);
 
743
  sc::auto_vec<double> xyz(new double[pc.size()*3]);
 
744
  for (int i=0; i<pc.size(); i++) {
 
745
    positions.get()[i] = &xyz.get()[i*3];
 
746
  }
 
747
 
 
748
  for (int ii=0; ii<pc.size(); ii++) {
 
749
    int i=pc[ii];
 
750
    charges.get()[ii] = molecule()->charge(i);
 
751
    for (int j=0; j<3; j++) positions.get()[ii][j] = molecule()->r(i,j);
 
752
  }
 
753
 
 
754
  Ref<PointChargeData> pcdata = new PointChargeData(pc.size(),
 
755
                                                    positions.get(),
 
756
                                                    charges.get());
 
757
  Ref<OneBodyOneCenterInt> ob = integral()->point_charge1(pcdata);
 
758
 
 
759
  const double *buffer = ob->buffer();
 
760
 
 
761
  for (int ii=0,icoef=0; ii<cd.size(); ii++) {
 
762
    int icenter = cd[ii];
 
763
    int joff = atom_basis()->shell_on_center(icenter,0);
 
764
    for (int j=0; j<atom_basis()->nshell_on_center(icenter); j++) {
 
765
      int jsh = j + joff;
 
766
      ob->compute_shell(jsh);
 
767
      int nfunc = atom_basis()->shell(jsh).nfunction();
 
768
      for (int k=0; k<nfunc; k++,icoef++) {
 
769
        e -= atom_basis_coef_[icoef] * buffer[k];
 
770
      }
 
771
    }
 
772
  }
 
773
 
 
774
  integral()->set_basis(basis());
 
775
 
 
776
  return e;
 
777
}
 
778
 
 
779
double
 
780
Wavefunction::nuc_rep_cd_cd(const std::vector<int>&c1,
 
781
                            const std::vector<int>&c2,
 
782
                            bool uniq)
 
783
{
 
784
  double e = 0.0;
 
785
 
 
786
  if (c1.size() == 0 || c2.size() == 0) return e;
 
787
 
 
788
  integral()->set_basis(atom_basis());
 
789
 
 
790
  Ref<TwoBodyTwoCenterInt> tb = integral()->electron_repulsion2();
 
791
 
 
792
  const double *buffer = tb->buffer();
 
793
 
 
794
  for (int ii=0; ii<c1.size(); ii++) {
 
795
    int icenter = c1[ii];
 
796
    int inshell = atom_basis()->nshell_on_center(icenter);
 
797
    int ish = atom_basis()->shell_on_center(icenter,0);
 
798
    for (int iish=0; iish<inshell; iish++,ish++) {
 
799
      int infunc = atom_basis()->shell(ish).nfunction();
 
800
      int ifuncoff = atom_basis()->shell_to_function(ish);
 
801
      int jjfence = (uniq?ii:c2.size());
 
802
      for (int jj=0; jj<jjfence; jj++) {
 
803
        int jcenter = c2[jj];
 
804
        int jnshell = atom_basis()->nshell_on_center(jcenter);
 
805
        int jsh = atom_basis()->shell_on_center(jcenter,0);
 
806
        for (int jjsh=0; jjsh<jnshell; jjsh++,jsh++) {
 
807
          int jnfunc = atom_basis()->shell(jsh).nfunction();
 
808
          int jfuncoff = atom_basis()->shell_to_function(jsh);
 
809
          tb->compute_shell(ish,jsh);
 
810
          for (int ifunc=0, ijfunc=0; ifunc<infunc; ifunc++) {
 
811
            for (int jfunc=0; jfunc<jnfunc; jfunc++, ijfunc++) {
 
812
              e += atom_basis_coef_[ifuncoff+ifunc]
 
813
                * atom_basis_coef_[jfuncoff+jfunc]
 
814
                * buffer[ijfunc];
 
815
            }
 
816
          }
 
817
        }
 
818
      }
 
819
    }
 
820
  }
 
821
 
 
822
  integral()->set_basis(basis());
 
823
 
 
824
  return e;
 
825
}
 
826
 
 
827
void
 
828
Wavefunction::nuclear_repulsion_energy_gradient(double *g)
 
829
{
 
830
  int natom = molecule()->natom();
 
831
  sc::auto_vec<double*> gp(new double*[natom]);
 
832
  for (int i=0; i<natom; i++) {
 
833
    gp.get()[i] = &g[i*3];
 
834
  }
 
835
  nuclear_repulsion_energy_gradient(gp.get());
 
836
}
 
837
 
 
838
void
 
839
Wavefunction::nuclear_repulsion_energy_gradient(double **g)
 
840
{
 
841
  if (atom_basis_.null()) {
 
842
    int natom = molecule()->natom();
 
843
    for (int i=0; i<natom; i++) {
 
844
      molecule()->nuclear_repulsion_1der(i,g[i]);
 
845
    }
 
846
    return;
 
847
  }
 
848
 
 
849
  // zero the gradient
 
850
  for (int i=0; i<molecule()->natom(); i++) {
 
851
    for (int j=0; j<3; j++) g[i][j] = 0.0;
 
852
  }
 
853
 
 
854
  // compute charge types
 
855
  std::vector<int> q_pc, q_cd, n_pc, n_cd;
 
856
  set_up_charge_types(q_pc,q_cd,n_pc,n_cd);
 
857
 
 
858
  // compute point charge - point charge contribution
 
859
  nuc_rep_grad_pc_pc(g, n_pc, n_pc, true /* i > j  */);
 
860
  nuc_rep_grad_pc_pc(g, q_pc, n_pc, false /* all i j */);
 
861
  if (molecule()->include_qq()) {
 
862
    nuc_rep_grad_pc_pc(g, q_pc, q_pc, true /* i > j */);
 
863
  }
 
864
 
 
865
  // compute point charge - charge distribution contribution
 
866
  nuc_rep_grad_pc_cd(g, n_pc, n_cd);
 
867
  nuc_rep_grad_pc_cd(g, q_pc, n_cd);
 
868
  nuc_rep_grad_pc_cd(g, n_pc, q_cd);
 
869
  if (molecule()->include_qq()) {
 
870
    nuc_rep_grad_pc_cd(g, q_pc, q_cd);
 
871
  }
 
872
 
 
873
  // compute the charge distribution - charge distribution contribution
 
874
  nuc_rep_grad_cd_cd(g, n_cd, n_cd, true /* i > j  */);
 
875
  nuc_rep_grad_cd_cd(g, q_cd, n_cd, false /* all i j  */);
 
876
  if (molecule()->include_qq()) {
 
877
    nuc_rep_grad_cd_cd(g, q_cd, q_cd, true /* i > j */);
 
878
  }
 
879
 
 
880
  // note: the electronic terms still need to be done in
 
881
  // a new hcore_deriv implemented in Wavefunction.
 
882
  throw std::runtime_error("Wavefunction::nuclear_repulsion_energy_gradient: not done");
 
883
 
 
884
}
 
885
 
 
886
void
 
887
Wavefunction::nuc_rep_grad_pc_pc(double **grad,
 
888
                                 const std::vector<int>&c1,
 
889
                                 const std::vector<int>&c2,
 
890
                                 bool uniq)
 
891
{
 
892
  if (c1.size() == 0 || c2.size() == 0) return;
 
893
 
 
894
  for (int ii=0; ii<c1.size(); ii++) {
 
895
    int i = c1[ii];
 
896
    SCVector3 ai(molecule()->r(i));
 
897
    double Zi = molecule()->charge(i);
 
898
    int jfence = (uniq?ii:c2.size());
 
899
    for (int jj=0; jj<jfence; jj++) {
 
900
      int j = c2[jj];
 
901
      SCVector3 aj(molecule()->r(j));
 
902
      double Zj = molecule()->charge(j);
 
903
      SCVector3 rdiff = ai - aj;
 
904
      double r2 = rdiff.dot(rdiff);
 
905
      double factor = - Zi * Zj / (r2*sqrt(r2));
 
906
      for (int k=0; k<3; k++) {
 
907
        grad[i][k] += factor * rdiff[k];
 
908
        grad[j][k] -= factor * rdiff[k];
 
909
      }
 
910
    }
 
911
  }
 
912
}
 
913
 
 
914
void
 
915
Wavefunction::nuc_rep_grad_pc_cd(double **grad,
 
916
                                 const std::vector<int>&pc,
 
917
                                 const std::vector<int>&cd)
 
918
{
 
919
  if (pc.size() == 0 || cd.size() == 0) return;
 
920
 
 
921
  throw std::runtime_error("Wavefunction::nuclear_repulsion_energy_gradient: not done");
 
922
}
 
923
 
 
924
void
 
925
Wavefunction::nuc_rep_grad_cd_cd(double **grad,
 
926
                                 const std::vector<int>&c1,
 
927
                                 const std::vector<int>&c2,
 
928
                                 bool uniq)
 
929
{
 
930
  if (c1.size() == 0 || c2.size() == 0) return;
 
931
 
 
932
  throw std::runtime_error("Wavefunction::nuclear_repulsion_energy_gradient: not done");
 
933
}
 
934
 
462
935
// returns the orthogonalization matrix
463
936
RefSCMatrix
464
937
Wavefunction::so_to_orthog_so()
480
953
  return gbs_;
481
954
}
482
955
 
 
956
Ref<Molecule>
 
957
Wavefunction::molecule() const
 
958
{
 
959
  return basis()->molecule();
 
960
}
 
961
 
 
962
Ref<GaussianBasisSet>
 
963
Wavefunction::atom_basis() const
 
964
{
 
965
  return atom_basis_;
 
966
}
 
967
 
 
968
const double *
 
969
Wavefunction::atom_basis_coef() const
 
970
{
 
971
  return atom_basis_coef_;
 
972
}
 
973
 
483
974
Ref<Integral>
484
975
Wavefunction::integral()
485
976
{
515
1006
Wavefunction::print(ostream&o) const
516
1007
{
517
1008
  MolecularEnergy::print(o);
 
1009
  ExEnv::out0() << indent << "Electronic basis:" << std::endl;
 
1010
  ExEnv::out0() << incindent;
518
1011
  basis()->print_brief(o);
 
1012
  ExEnv::out0() << decindent;
 
1013
  if (atom_basis_.nonnull()) {
 
1014
    ExEnv::out0() << indent << "Nuclear basis:" << std::endl;
 
1015
    ExEnv::out0() << incindent;
 
1016
    atom_basis_->print_brief(o);
 
1017
    ExEnv::out0() << decindent;
 
1018
  }
519
1019
  // the other stuff is a wee bit too big to print
520
1020
  if (print_nao_ || print_npa_) {
521
1021
    tim_enter("NAO");
614
1114
  return orthog_method_;
615
1115
}
616
1116
 
 
1117
void
 
1118
Wavefunction::set_orthog_method(const OverlapOrthog::OrthogMethod& omethod)
 
1119
{
 
1120
  if (orthog_method_ != omethod) {
 
1121
    orthog_method_ = omethod;
 
1122
    init_orthog();
 
1123
    obsolete();
 
1124
  }
 
1125
}
 
1126
 
617
1127
double
618
1128
Wavefunction::lindep_tol() const
619
1129
{
620
1130
  return lindep_tol_;
621
1131
}
622
1132
 
 
1133
void
 
1134
Wavefunction::set_lindep_tol(double lindep_tol)
 
1135
{
 
1136
  if (lindep_tol_ != lindep_tol) {
 
1137
    lindep_tol_ = lindep_tol;
 
1138
    init_orthog();
 
1139
    obsolete();
 
1140
  }
 
1141
}
 
1142
 
 
1143
void
 
1144
Wavefunction::scale_atom_basis_coef()
 
1145
{
 
1146
  std::vector<int> i_cd;
 
1147
  for (int i=0; i<atom_basis_->ncenter(); i++) {
 
1148
    if (atom_basis_->nshell_on_center(i) > 0) i_cd.push_back(i);
 
1149
  }
 
1150
 
 
1151
  if (atom_basis_->max_angular_momentum() > 0) {
 
1152
    // Only s functions will preserve the full symmetry.
 
1153
    // Can only normalize functions that don't integrate to zero.
 
1154
    throw std::runtime_error("ChargeDistInt: max am larger than 0");
 
1155
  }
 
1156
 
 
1157
  int coef_offset = 0;
 
1158
  int icoef = 0;
 
1159
  for (int ii=0; ii<i_cd.size(); ii++) {
 
1160
    int i = i_cd[ii];
 
1161
    int nshell = atom_basis_->nshell_on_center(i);
 
1162
    int nfunc = 0;
 
1163
    if (nshell > 0) {
 
1164
      double raw_charge = 0.0;
 
1165
      for (int jj=0, icoef=coef_offset; jj<nshell; jj++) {
 
1166
        int j = atom_basis_->shell_on_center(i,jj);
 
1167
        const GaussianShell &shell = atom_basis_->shell(j);
 
1168
        // loop over the contractions
 
1169
        // the number of functions in each contraction is 1
 
1170
        // since only s functions are allowed
 
1171
        for (int k=0; k<shell.ncontraction(); k++, icoef++) {
 
1172
          for (int l=0; l<shell.nprimitive(); l++) {
 
1173
            double alpha = shell.exponent(l);
 
1174
            double con_coef = shell.coefficient_unnorm(k,l);
 
1175
            // The exponent is halved because the normalization
 
1176
            // formula is for the s function squared.
 
1177
            double integral = ::pow(M_PI/alpha,1.5);
 
1178
            raw_charge += atom_basis_coef_[icoef] * con_coef * integral;
 
1179
//             std::cout << "con_coef = " << con_coef
 
1180
//                       << " integral = " << integral
 
1181
//                       << std::endl;
 
1182
          }
 
1183
        }
 
1184
        nfunc += shell.ncontraction();
 
1185
      }
 
1186
      double charge = atom_basis_->molecule()->charge(i);
 
1187
      double correction = charge/raw_charge;
 
1188
      for (int icoef=coef_offset; icoef<coef_offset+nfunc; icoef++) {
 
1189
        atom_basis_coef_[icoef] = correction * atom_basis_coef_[icoef];
 
1190
      }
 
1191
    }
 
1192
    coef_offset += nshell;
 
1193
  }
 
1194
}
 
1195
 
623
1196
/////////////////////////////////////////////////////////////////////////////
624
1197
 
625
1198
// Local Variables: