btest.cc

大型并行量子化学软件；支持密度泛函（DFT）。可以进行各种量子化学计算。支持CHARMM并行计算。非常具有应用价值。

  double *b_val_plszx = new double[nbasis];
  double *b_val_mnszx = new double[nbasis];
  double *g_val = new double[3*nbasis];
  double *h_val = new double[6*nbasis];

  const int x_ = 0;
  const int y_ = 1;
  const int z_ = 2;
  const int xx_ = 0;
  const int yx_ = 1;
  const int yy_ = 2;
  const int zx_ = 3;
  const int zy_ = 4;
  const int zz_ = 5;

  SCVector3 r;
  SCVector3 d;
  double delta = 0.001;
  SCVector3 dx(delta, 0., 0.);
  SCVector3 dy(0., delta, 0.);
  SCVector3 dz(0., 0., delta);
  SCVector3 dxy(delta, delta, 0.);
  SCVector3 dxz(delta, 0., delta);
  SCVector3 dyz(0., delta, delta);
  double deltax = 0.1;
  GaussianBasisSet::ValueData vdat(gbs, integral);
  for (r.x()=0.0; r.x() < 1.0; r.x() += deltax) {
      deltax *= 2.;
      double deltay = 0.1;
      for (r.y()=0.0; r.y() < 1.0; r.y() += deltay) {
          deltay *= 2.;
          double deltaz = 0.1;
          for (r.z()=0.0; r.z() < 1.0; r.z() += deltaz) {
              deltaz *= 2.;
              cout << "R = " << r << endl;
              gbs->hessian_values(r, &vdat, h_val, g_val, b_val);
              gbs->values(r + dx, &vdat, b_val_plsx);
              gbs->values(r - dx, &vdat, b_val_mnsx);
              gbs->values(r + dy, &vdat, b_val_plsy);
              gbs->values(r - dy, &vdat, b_val_mnsy);
              gbs->values(r + dz, &vdat, b_val_plsz);
              gbs->values(r - dz, &vdat, b_val_mnsz);
              gbs->values(r + dxy, &vdat, b_val_plsyx);
              gbs->values(r - dxy, &vdat, b_val_mnsyx);
              gbs->values(r + dyz, &vdat, b_val_plszy);
              gbs->values(r - dyz, &vdat, b_val_mnszy);
              gbs->values(r + dxz, &vdat, b_val_plszx);
              gbs->values(r - dxz, &vdat, b_val_mnszx);
              for (int i=0; i<nbasis; i++) {
                  int shell = gbs->function_to_shell(i);
                  int func = i - gbs->shell_to_function(shell);
                  double g_val_test[3];
                  double h_val_test[6];
                  int x = i*3+x_;
                  int y = i*3+y_;
                  int z = i*3+z_;
                  g_val_test[x_] = 0.5*(b_val_plsx[i] - b_val_mnsx[i])/delta;
                  g_val_test[y_] = 0.5*(b_val_plsy[i] - b_val_mnsy[i])/delta;
                  g_val_test[z_] = 0.5*(b_val_plsz[i] - b_val_mnsz[i])/delta;
                  int xx = i*6+xx_;
                  int yx = i*6+yx_;
                  int yy = i*6+yy_;
                  int zx = i*6+zx_;
                  int zy = i*6+zy_;
                  int zz = i*6+zz_;
                  h_val_test[xx_]
                      = (b_val_plsx[i] + b_val_mnsx[i] - 2. * b_val[i])
                        * 1./(delta*delta);
                  h_val_test[yy_]
                      = (b_val_plsy[i] + b_val_mnsy[i] - 2. * b_val[i])
                        * 1./(delta*delta);
                  h_val_test[zz_]
                      = (b_val_plsz[i] + b_val_mnsz[i] - 2. * b_val[i])
                        * 1./(delta*delta);
                  h_val_test[yx_]
                      = 0.5 * ((b_val_plsyx[i]+b_val_mnsyx[i]-2.*b_val[i])
                               * 1. / (delta * delta)
                               - h_val_test[xx_] - h_val_test[yy_]);
                  h_val_test[zx_]
                      = 0.5 * ((b_val_plszx[i]+b_val_mnszx[i]-2.*b_val[i])
                               * 1. / (delta * delta)
                               - h_val_test[xx_] - h_val_test[zz_]);
                  h_val_test[zy_]
                      = 0.5 * ((b_val_plszy[i]+b_val_mnszy[i]-2.*b_val[i])
                               * 1. / (delta * delta)
                               - h_val_test[zz_] - h_val_test[yy_]);
                  checkerror("x", shell, func, g_val_test[x_], g_val[x]);
                  checkerror("y", shell, func, g_val_test[y_], g_val[y]);
                  checkerror("z", shell, func, g_val_test[z_], g_val[z]);
                  checkerror("xx", shell, func, h_val_test[xx_], h_val[xx]);
                  checkerror("yy", shell, func, h_val_test[yy_], h_val[yy]);
                  checkerror("zz", shell, func, h_val_test[zz_], h_val[zz]);
                  checkerror("yx", shell, func, h_val_test[yx_], h_val[yx]);
                  checkerror("zx", shell, func, h_val_test[zx_], h_val[zx]);
                  checkerror("zy", shell, func, h_val_test[zy_], h_val[zy]);
                }
            }
        }
    }
  delete[] b_val;
  delete[] b_val_plsx;
  delete[] b_val_mnsx;
  delete[] b_val_plsy;
  delete[] b_val_mnsy;
  delete[] b_val_plsz;
  delete[] b_val_mnsz;
  delete[] b_val_plsyx;
  delete[] b_val_mnsyx;
  delete[] b_val_plszy;
  delete[] b_val_mnszy;
  delete[] b_val_plszx;
  delete[] b_val_mnszx;
  delete[] g_val;
  delete[] h_val;
}

void
do_extent_test(const Ref<GaussianBasisSet> &gbs)
{
  int i, j;
  for (i=0; i<gbs->nshell(); i++) {
      gbs->shell(i).print();
      for (j=0; j<10; j++) {
          cout << " " << gbs->shell(i).monobound(0.1*j);
        }
      cout << endl;
      for (j=0; j<10; j++) {
          double threshold = pow(10.0, -j);
          //cout << " threshold = " << threshold << endl;
          cout << " " << gbs->shell(i).extent(threshold);
        }
      cout << endl;
    }

  Ref<ShellExtent> extent = new ShellExtent;
  extent->init(gbs);
  extent->print();
}

void
do_gpetite_test(const Ref<GaussianBasisSet> &b1,
                const Ref<GaussianBasisSet> &b2)
{
  canonical_aabc c4(b1,b1,b1,b2);
  Ref<GPetite4<canonical_aabc> > p4
      = new GPetite4<canonical_aabc>(b1,b1,b1,b2,c4);
  cout << "tesing GPetite4<canonical_aabc>" << endl;
  for (int i=0; i<b1->nshell(); i++) {
      for (int j=0; j<=i; j++) {
          for (int k=0; k<b1->nshell(); k++) {
              for (int l=0; l<b2->nshell(); l++) {
                  cout << " " << i
                       << " " << j
                       << " " << k
                       << " " << l
                       << " in p4: " << p4->in_p4(i,j,k,l)
                       << endl;
                }
            }
        }
    }
}

int
main(int, char *argv[])
{
  int i, j;

  char o[10000];
#ifdef HAVE_SSTREAM
  ostringstream perlout(o);
#else
  ostrstream perlout(o,sizeof(o));
#endif

  const char *filename = (argv[1]) ? argv[1] : SRCDIR "/btest.kv";
  
  Ref<KeyVal> keyval = new ParsedKeyVal(filename);
  
  Ref<Integral> intgrl = new IntegralV3;

  int dooverlap = keyval->booleanvalue("overlap");
  int doeigvals = keyval->booleanvalue("eigvals");
  int dostate = keyval->booleanvalue("state");
  int doso = keyval->booleanvalue("so");
  int doatoms = keyval->booleanvalue("atoms");
  int dopetite = keyval->booleanvalue("petite");
  int dogpetite = keyval->booleanvalue("gpetite");
  int dovalues = keyval->booleanvalue("values");
  int doextent = keyval->booleanvalue("extent");
  int doaoorthog = keyval->booleanvalue("aoorthog");

  for (i=0; i<keyval->count("test"); i++) {
      Ref<GaussianBasisSet> gbs;
      gbs << keyval->describedclassvalue("test", i);
      Ref<GaussianBasisSet> gbs2;
      gbs2 << keyval->describedclassvalue("test2", i);

      if (dooverlap) test_overlap(gbs,gbs2,intgrl);

      if (doaoorthog) test_aoorthog(gbs,intgrl);

      if (doeigvals) test_eigvals(gbs,intgrl);

      if (dostate) {
          StateOutBin out("btest.out");
          SavableState::save_state(gbs.pointer(),out);
          out.close();
          StateInBin in("btest.out");
          gbs << SavableState::restore_state(in);
          gbs->print();
        }

      if (dopetite) {
          intgrl->set_basis(gbs);
          intgrl->petite_list()->print();
        }

      if (dogpetite) {
          do_gpetite_test(gbs,gbs2);
        }

      if (doso) {
          do_so_test(keyval, intgrl, gbs);
        }

      if (dovalues) {
          intgrl->set_basis(gbs);
          test_func_values(gbs,intgrl);
        }

      if (doextent) {
          do_extent_test(gbs);
        }
    }

  if (doatoms) {
      const int nelem = 37;

      // Make H, C, and P molecules
      Ref<Molecule> hmol = new Molecule();
      hmol->add_atom(AtomInfo::string_to_Z("H"),0,0,0);
      Ref<Molecule> cmol = new Molecule();
      cmol->add_atom(AtomInfo::string_to_Z("C"),0,0,0);
      Ref<Molecule> pmol = new Molecule();
      pmol->add_atom(AtomInfo::string_to_Z("P"),0,0,0);

      perlout << "%basissets = (" << endl;
      int nbasis = keyval->count("basislist");
      Ref<KeyVal> nullkv = new AssignedKeyVal();
      for (i=0; i<nbasis; i++) {
          int first_element = 1;
          char *basisname = keyval->pcharvalue("basislist",i);
          perlout << "  \"" << basisname << "\" => (";
          BasisFileSet bfs(nullkv);
          Ref<KeyVal> basiskv = bfs.keyval(nullkv, basisname);
          char elemstr[512];
          elemstr[0] = '\0';
          int last_elem_exists = 0;
          int n0 = 0;
          int n1 = 0;
          int n2 = 0;
          for (j=0; j<nelem; j++) {
              Ref<AssignedKeyVal> atombaskv_a(new AssignedKeyVal());
              Ref<KeyVal> atombaskv(atombaskv_a);
              char keyword[256];
              strcpy(keyword,":basis:");
              strcat(keyword,AtomInfo::name(j+1));
              strcat(keyword,":");
              strcat(keyword,basisname);
              if (basiskv->exists(keyword)) {
                  if (!first_element) {
                      perlout << ",";
                    }
                  else {
                      first_element = 0;
                    }
                  perlout << "\"" << AtomInfo::symbol(j+1) << "\"";
                  if (!last_elem_exists) {
                      if (elemstr[0] != '\0') strcat(elemstr,", ");
                      strcat(elemstr,AtomInfo::symbol(j+1));
                    }
                  else if (last_elem_exists == 2) {
                      strcat(elemstr,"-");
                    }
                  last_elem_exists++;
                  if (j+1 == 1) {
                      atombaskv_a->assign("name", basisname);
                      atombaskv_a->assign("molecule", hmol.pointer());
                      Ref<GaussianBasisSet> gbs=new GaussianBasisSet(atombaskv);
                      n0 = gbs->nbasis();
                    }
                  if (j+1 == 6) {
                      atombaskv_a->assign("name", basisname);
                      atombaskv_a->assign("molecule", cmol.pointer());
                      Ref<GaussianBasisSet> gbs=new GaussianBasisSet(atombaskv);
                      n1 = gbs->nbasis();
                    }
                  if (j+1 == 15) {
                      atombaskv_a->assign("name", basisname);
                      atombaskv_a->assign("molecule", pmol.pointer());
                      Ref<GaussianBasisSet> gbs=new GaussianBasisSet(atombaskv);
                      n2 = gbs->nbasis();
                    }
                }
              else {
                  if (last_elem_exists > 1) {
                      if (last_elem_exists == 2) strcat(elemstr,", ");
                      strcat(elemstr, AtomInfo::symbol(j));
                    }
                  last_elem_exists = 0;
                }
            }
          perlout << ")";
          if (i != nbasis-1) perlout << "," << endl;
          perlout << endl;
          cout << "<tr><td><tt>" << basisname
               << "</tt><td>" << elemstr << "<td>";
          if (n0>0) cout << n0;
          cout << "<td>";
          if (n1>0) cout << n1;
          cout << "<td>";
          if (n2>0) cout << n2;
          cout << endl;
          delete[] basisname;
        }

      perlout << ")" << endl;

#ifdef HAVE_SSTREAM
      const char *perlout_s = perlout.str().c_str();
#else
      perlout << ")" << ends;
      char *perlout_s = perlout.str();
#endif
      cout << perlout_s;
    }

  return 0;
}

/////////////////////////////////////////////////////////////////////////////

// Local Variables:
// mode: c++
// c-file-style: "CLJ"
// End: