Add functionality to compute energy and forces

src/DensityMatrix.cc

@@ -149,14 +149,11 @@ template <class MatrixType>
 void DensityMatrix<MatrixType>::setUniform(
     const double nel, const int new_orbitals_index)
 {
-#ifdef PRINT_OPERATIONS
-    MGmol_MPI& mmpi = *(MGmol_MPI::instance());
+    MGmol_MPI& mmpi  = *(MGmol_MPI::instance());
+    const double occ = (double)((double)nel / (double)dim_);
     if (mmpi.instancePE0())
-        std::cout << "template <class MatrixType> "
-                     "DensityMatrix<MatrixType>::setUniform()"
+        std::cout << "DensityMatrix::setUniform(), occupation = " << occ
                   << std::endl;
-#endif
-    const double occ = (double)((double)nel / (double)dim_);
     assert(occ < 1.01);
     for (int i = 0; i < dim_; i++)
         occupation_[i] = occ;

src/MGmol.cc

@@ -1432,15 +1432,15 @@ void MGmol<OrbitalsType>::getAtomicNumbers(std::vector<short>& an)
 
 template <class OrbitalsType>
 double MGmol<OrbitalsType>::evaluateEnergyAndForces(
-    const std::vector<double>& tau, std::vector<short>& atnumbers,
+    const std::vector<double>& tau, const std::vector<short>& atnumbers,
     std::vector<double>& forces)
 {
     return evaluateEnergyAndForces(current_orbitals_, tau, atnumbers, forces);
 }
 
 template <class OrbitalsType>
 double MGmol<OrbitalsType>::evaluateEnergyAndForces(Orbitals* orbitals,
-    const std::vector<double>& tau, std::vector<short>& atnumbers,
+    const std::vector<double>& tau, const std::vector<short>& atnumbers,
     std::vector<double>& forces)
 {
     assert(tau.size() == 3 * atnumbers.size());
@@ -1462,6 +1462,46 @@ double MGmol<OrbitalsType>::evaluateEnergyAndForces(Orbitals* orbitals,
     return eks;
 }
 
+template <class OrbitalsType>
+double MGmol<OrbitalsType>::evaluateDMandEnergyAndForces(Orbitals* orbitals,
+    const std::vector<double>& tau, const std::vector<short>& atnumbers,
+    std::vector<double>& forces)
+{
+    OrbitalsType* dorbitals = dynamic_cast<OrbitalsType*>(orbitals);
+
+    ions_->setPositions(tau, atnumbers);
+
+    moveVnuc(*ions_);
+
+    // initialize electronic density
+    rho_->update(*dorbitals);
+
+    // initialize potential
+    update_pot(*ions_);
+
+    // initialize projected matrices
+    updateHmatrix(*dorbitals, *ions_);
+    proj_matrices_->updateThetaAndHB();
+
+    // compute DM
+    std::shared_ptr<DMStrategy<OrbitalsType>> dm_strategy(
+        DMStrategyFactory<OrbitalsType,
+            dist_matrix::DistMatrix<double>>::create(comm_, os_, *ions_,
+            rho_.get(), energy_.get(), electrostat_.get(), this,
+            proj_matrices_.get(), dorbitals));
+
+    dm_strategy->update(*dorbitals);
+
+    // evaluate energy and forces
+    double ts = 0.;
+    double eks
+        = energy_->evaluateTotal(ts, proj_matrices_.get(), *dorbitals, 2, os_);
+
+    force(*dorbitals, *ions_);
+
+    return eks;
+}
+
 template class MGmol<LocGridOrbitals>;
 template class MGmol<ExtendedGridOrbitals>;
 template int MGmol<LocGridOrbitals>::initial<MemorySpace::Host>();

src/MGmol.h

@@ -186,11 +186,30 @@ class MGmol : public MGmolInterface
 
     void run() override;
 
+    /*
+     * Evaluate the energy and forces for an atomic configuration
+     * specified by tau (input)
+     */
     double evaluateEnergyAndForces(const std::vector<double>& tau,
-        std::vector<short>& atnumbers, std::vector<double>& forces);
+        const std::vector<short>& atnumbers, std::vector<double>& forces);
 
-    double evaluateEnergyAndForces(Orbitals*, const std::vector<double>& tau,
-        std::vector<short>& atnumbers, std::vector<double>& forces);
+    /*
+     * Evaluate the energy and forces for an atomic configuration
+     * specified by tau (input), using the input "orbitals" as initial
+     * guess for the wavefunctions
+     */
+    double evaluateEnergyAndForces(Orbitals* orbitals,
+        const std::vector<double>& tau, const std::vector<short>& atnumbers,
+        std::vector<double>& forces);
+
+    /*
+     * Evaluate the energy and forces for an atomic configuration
+     * specified by tau (input), using the input "orbitals" as wavefunctions
+     * (fixed)
+     */
+    double evaluateDMandEnergyAndForces(Orbitals* orbitals,
+        const std::vector<double>& tau, const std::vector<short>& atnumbers,
+        std::vector<double>& forces);
 
     /*
      * get internal atomic positions

src/MGmolInterface.h

@@ -25,13 +25,19 @@ class MGmolInterface
     virtual int setupConstraintsFromInput(const std::string input_file) = 0;
     virtual void setup()                                                = 0;
     virtual void run()                                                  = 0;
+
     virtual double evaluateEnergyAndForces(const std::vector<double>& tau,
-        std::vector<short>& atnumbers, std::vector<double>& forces)
+        const std::vector<short>& atnumbers, std::vector<double>& forces)
         = 0;
     virtual double evaluateEnergyAndForces(Orbitals*,
-        const std::vector<double>& tau, std::vector<short>& atnumbers,
+        const std::vector<double>& tau, const std::vector<short>& atnumbers,
+        std::vector<double>& forces)
+        = 0;
+    virtual double evaluateDMandEnergyAndForces(Orbitals*,
+        const std::vector<double>& tau, const std::vector<short>& atnumbers,
         std::vector<double>& forces)
         = 0;
+
     virtual void getAtomicPositions(std::vector<double>& tau) = 0;
     virtual void getAtomicNumbers(std::vector<short>& an)     = 0;
     virtual std::shared_ptr<ProjectedMatricesInterface> getProjectedMatrices()

src/MVPSolver.cc

@@ -208,6 +208,7 @@ int MVPSolver<OrbitalsType, MatrixType>::solve(OrbitalsType& orbitals)
 
         kbpsi.computeHvnlMatrix(&kbpsi, ions_, h11_nl);
 
+        const double tol_de0 = 1.e-12;
         for (int inner_it = 0; inner_it < n_inner_steps_; inner_it++)
         {
             if (onpe0 && ct.verbose > 1)
@@ -270,6 +271,14 @@ int MVPSolver<OrbitalsType, MatrixType>::solve(OrbitalsType& orbitals)
 
                 double de0 = evaluateDerivative(dmInit, delta_dm, ts0);
 
+                if (std::abs(de0) < tol_de0 && inner_it > 0)
+                {
+                    if (onpe0 && ct.verbose > 0)
+                        std::cout << "MVP: de0 = " << de0
+                                  << ", convergence achieved" << std::endl;
+                    break;
+                }
+
                 //
                 // evaluate free energy at beta=1
                 //
@@ -306,6 +315,7 @@ int MVPSolver<OrbitalsType, MatrixType>::solve(OrbitalsType& orbitals)
                 // line minimization
                 const double beta
                     = minQuadPolynomial(e0, e1, de0, (ct.verbose > 2), os_);
+                assert(!std::isnan(beta));
 
                 if (onpe0 && ct.verbose > 0)
                 {
@@ -348,6 +358,8 @@ int MVPSolver<OrbitalsType, MatrixType>::solve(OrbitalsType& orbitals)
                 orbitals.getProjMatrices());
         projmatrices->setDM(*work_, orbitals.getIterativeIndex());
         projmatrices->setEigenvalues(proj_mat_work_->getEigenvalues());
+        projmatrices->assignH(proj_mat_work_->getH());
+        projmatrices->setHB2H();
     }
 
     // Generate new density

src/ProjectedMatrices.cc

@@ -650,8 +650,9 @@ double ProjectedMatrices<MatrixType>::computeEntropy()
         else
         {
             if (mmpi.PE0() && ct.verbose > 1)
-                (*MPIdata::sout)
-                    << "occupations uptodate, skip computation..." << std::endl;
+                (*MPIdata::sout) << "computeEntropy: occupations uptodate, "
+                                    "skip computation..."
+                                 << std::endl;
         }
         entropy = computeEntropy(width_);
     }

src/Rho.cc

@@ -90,16 +90,16 @@ void Rho<OrbitalsType>::update(OrbitalsType& current_orbitals)
 
     update_tm_.start();
 
-    if (verbosity_level_ > 2 && onpe0)
-        (*MPIdata::sout) << "Rho<OrbitalsType>::update()" << std::endl;
+    if (verbosity_level_ > 1 && onpe0)
+        (*MPIdata::sout) << "Rho::update()..." << std::endl;
 
     const int new_iterative_index
         = ((1 + current_orbitals.getIterativeIndex()) % 100)
           + (proj_matrices.getDMMatrixIndex() % 100) * 100;
 
     if (iterative_index_ == new_iterative_index)
     {
-        if (onpe0 && verbosity_level_ > 2)
+        if (onpe0 && verbosity_level_ > 1)
             (*MPIdata::sout) << "Rho already up to date, iterative_index_="
                              << iterative_index_ << std::endl;
         return;

tests/CMakeLists.txt

@@ -250,6 +250,8 @@ add_executable(testEnergyAndForces
                ${CMAKE_SOURCE_DIR}/tests/EnergyAndForces/testEnergyAndForces.cc)
 add_executable(testWFEnergyAndForces
                ${CMAKE_SOURCE_DIR}/tests/WFEnergyAndForces/testWFEnergyAndForces.cc)
+add_executable(testDMandEnergyAndForces
+               ${CMAKE_SOURCE_DIR}/tests/DMandEnergyAndForces/testDMandEnergyAndForces.cc)
 
 if(${MAGMA_FOUND})
   add_executable(testOpenmpOffload
@@ -353,6 +355,14 @@ add_test(NAME testWFEnergyAndForces
          ${CMAKE_CURRENT_SOURCE_DIR}/WFEnergyAndForces/mgmol.cfg
          ${CMAKE_CURRENT_SOURCE_DIR}/WFEnergyAndForces/sih4.xyz
          ${CMAKE_CURRENT_SOURCE_DIR}/../potentials)
+add_test(NAME testDMandEnergyAndForces
+         COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/DMandEnergyAndForces/test.py
+         ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS}
+         ${CMAKE_CURRENT_BINARY_DIR}/testDMandEnergyAndForces
+         ${CMAKE_CURRENT_SOURCE_DIR}/DMandEnergyAndForces/mgmol.cfg
+         ${CMAKE_CURRENT_SOURCE_DIR}/DMandEnergyAndForces/coords.in
+         ${CMAKE_CURRENT_SOURCE_DIR}/DMandEnergyAndForces/lrs.in
+         ${CMAKE_CURRENT_SOURCE_DIR}/../potentials)
 
 if(${MAGMA_FOUND}) 
   add_test(NAME testOpenmpOffload
@@ -539,6 +549,7 @@ target_link_libraries(testSuperSampling PRIVATE MPI::MPI_CXX)
 target_link_libraries(testDirectionalReduce PRIVATE MPI::MPI_CXX)
 target_link_libraries(testEnergyAndForces PRIVATE mgmol_src)
 target_link_libraries(testWFEnergyAndForces PRIVATE mgmol_src)
+target_link_libraries(testDMandEnergyAndForces PRIVATE mgmol_src)
 
 if(${MAGMA_FOUND})
   target_link_libraries(testDistVector PRIVATE ${SCALAPACK_LIBRARIES}

tests/DMandEnergyAndForces/coords.in

@@ -0,0 +1,2 @@
+N1  1  0.  0.  -1.0345
+N2  1  0.  0.   1.0345

tests/DMandEnergyAndForces/lrs.in

@@ -0,0 +1,6 @@
+0.00    0.7    1.4
+0.00   -0.7   -1.4
+0.00    0.7   -1.4
+0.00   -0.7    1.4
+0.46    0.0    0.0
+-0.46    0.0    0.0

tests/DMandEnergyAndForces/mgmol.cfg

@@ -0,0 +1,35 @@
+verbosity=2
+xcFunctional=LDA
+FDtype=4th
+[Mesh]
+nx=64
+ny=64
+nz=64
+[Domain]
+ox=-6.
+oy=-6.
+oz=-6.
+lx=12.
+ly=12.
+lz=12.
+[Potentials]
+pseudopotential=pseudo.N_ONCVPSP_LDA
+[Run]
+type=QUENCH
+[Quench]
+solver=PSD
+max_steps=100
+atol=1.e-8
+step_length=2.
+ortho_freq=10
+[Orbitals]
+initial_type=Gaussian
+initial_width=1.5
+temperature=10.
+nempty=2
+[Restart]
+output_level=3
+output_filename=WF
+[DensityMatrix]
+solver=MVP
+nb_inner_it=1

tests/DMandEnergyAndForces/test.py

@@ -0,0 +1,73 @@
+#!/usr/bin/env python
+import sys
+import os
+import subprocess
+import string
+
+print("Test DMandEnergyAndForces...")
+
+nargs=len(sys.argv)
+
+mpicmd = sys.argv[1]+" "+sys.argv[2]+" "+sys.argv[3]
+for i in range(4,nargs-6):
+  mpicmd = mpicmd + " "+sys.argv[i]
+print("MPI run command: {}".format(mpicmd))
+
+exe = sys.argv[nargs-5]
+inp = sys.argv[nargs-4]
+coords = sys.argv[nargs-3]
+print("coordinates file: %s"%coords)
+lrs = sys.argv[-2]
+
+#create links to potentials files
+dst = 'pseudo.N_ONCVPSP_LDA'
+src = sys.argv[-1] + '/' + dst
+
+if not os.path.exists(dst):
+  print("Create link to %s"%dst)
+  os.symlink(src, dst)
+
+#run
+command = "{} {} -c {} -i {} -l {}".format(mpicmd,exe,inp,coords,lrs)
+print("Run command: {}".format(command))
+output = subprocess.check_output(command,shell=True)
+lines=output.split(b'\n')
+
+#analyse output
+energies=[]
+for line in lines:
+  if line.count(b'%%'):
+    print(line)
+    words=line.split()
+    words=words[5].split(b',')[0]
+    energy = words.decode()
+  if line.count(b'achieved'):
+    energies.append(energy)
+    break
+
+for line in lines:
+  if line.count(b'MVP') and line.count(b'iteration'):
+    print(line)
+  if line.count(b'Eks2'):
+    print(line)
+    words=line.split()
+    word=words[2]
+    energy = word.decode()
+    energies.append(energy)
+    break
+
+print("Check energies...")
+print( energies )
+if len(energies)<2:
+  print("Expected two converged energies")
+  sys.exit(1)
+
+tol = 1.e-6
+diff=eval(energies[1])-eval(energies[0])
+print(diff)
+if abs(diff)>tol:
+  print("Energies differ: {} vs {} !!!".format(energies[0],energies[1]))
+  sys.exit(1)
+
+print("Test SUCCESSFUL!")
+sys.exit(0)