Print gradient (#22)

src/program_dftd.cpp

@@ -24,32 +24,33 @@
 #include "dftd_damping.h"
 #include "dftd_dispersion.h"
 #include "dftd_geometry.h"
-#include "dftd_model.h"
 #include "dftd_matrix.h"
+#include "dftd_model.h"
 #include "dftd_readxyz.h"
 
 class argparser {
-    public:
-  argparser(int &argc, char **argv) {
-    for (int i = 1; i != argc; i++) {
-      this->args.push_back(std::string(argv[i]));
+  public:
+    argparser(int &argc, char **argv) {
+      for (int i = 1; i != argc; i++) {
+        this->args.push_back(std::string(argv[i]));
+      }
     }
-  }
-  const std::string &getopt(const std::string &opt) const {
-    std::vector<std::string>::const_iterator iter;
-    iter = find(this->args.begin(), this->args.end(), opt);
-    if (iter != this->args.end() && ++iter != this->args.end()) {
-      return *iter;
+    const std::string &getopt(const std::string &opt) const {
+      std::vector<std::string>::const_iterator iter;
+      iter = find(this->args.begin(), this->args.end(), opt);
+      if (iter != this->args.end() && ++iter != this->args.end()) {
+        return *iter;
+      }
+      static const std::string empty("");
+      return empty;
+    }
+    bool getflag(const std::string &opt) const {
+      return find(this->args.begin(), this->args.end(), opt) !=
+             this->args.end();
     }
-    static const std::string empty("");
-    return empty;
-  }
-  bool getflag(const std::string &opt) const {
-    return find(this->args.begin(), this->args.end(), opt) != this->args.end();
-  }
 
-    private:
-  std::vector<std::string> args;
+  private:
+    std::vector<std::string> args;
 };
 
 void dftd4_citation() {
@@ -125,13 +126,17 @@ int main(int argc, char **argv) {
     help();
     exit(EXIT_FAILURE);
   }
+
   // setup the argparser from the commandline
   argparser args(argc, argv);
+
   // check for help flag first
   if (args.getflag("-h") || args.getflag("--help")) {
     help();
     exit(EXIT_SUCCESS);
   }
+
+  // get other flags
   if (args.getflag("-v") || args.getflag("--verbose")) { lverbose = true; }
   if (args.getflag("-g") || args.getflag("--grad")) { lgrad = true; }
   if (args.getflag("--func")) {
@@ -149,10 +154,46 @@ int main(int argc, char **argv) {
   dftd4::TCutoff cutoff;
   dftd4::TD4Model d4;
 
-  info = dftd4::get_dispersion(mol, charge, d4, par, cutoff, energy, nullptr);
-  if (info != 0) return EXIT_FAILURE;
+  // masking (nothing excluded)
+  dftd4::TVector<int> realIdx;
+  realIdx.NewVec(mol.NAtoms);
+  int nat = 0;
+  for (int i = 0; i != mol.NAtoms; i++) {
+    realIdx(i) = nat;
+    nat++;
+  }
 
-  std::cout << "Dispersion energy: " << energy << " Eh\n";
+  // analytical gradient
+  double *d4grad;
+  if (lgrad) {
+    d4grad = new double[3 * mol.NAtoms];
+    for (int i = 0; i < 3 * mol.NAtoms; i++) {
+      d4grad[i] = 0.0;
+    }
+  } else {
+    d4grad = nullptr;
+  }
+
+  info = dftd4::get_dispersion(mol, charge, d4, par, cutoff, energy, d4grad);
+  if (info != EXIT_SUCCESS) return info;
+
+  // Print results
+  printf("Dispersion energy: %.15e Eh\n", energy);
+
+  if (lgrad) {
+    printf("\nDispersion gradient\n");
+    int count{0};
+    for (int i = 0; i < mol.NAtoms; i++) {
+      printf(
+        "%3i  :  %+14.9le  %+14.9le  %+14.9le\n",
+        i + 1,
+        d4grad[count],
+        d4grad[count + 1],
+        d4grad[count + 2]
+      );
+      count = count + 3;
+    };
+  }
 
   mol.FreeMemory();