Fix and test EnergyAndForces interface

src/MGmol.cc

@@ -1428,6 +1428,8 @@ double MGmol<OrbitalsType>::evaluateEnergyAndForces(
 
     ions_->setPositions(tau, atnumbers);
 
+    moveVnuc(*ions_);
+
     double eks = 0.;
     quench(current_orbitals_, *ions_, ct.max_electronic_steps, 20, eks);
 

src/MGmol.h

@@ -180,7 +180,10 @@ class MGmol : public MGmolInterface
 
     /* access functions */
     OrbitalsType* getOrbitals() { return current_orbitals_; }
-    std::shared_ptr<Hamiltonian<OrbitalsType>> getHamiltonian() { return hamiltonian_; }
+    std::shared_ptr<Hamiltonian<OrbitalsType>> getHamiltonian()
+    {
+        return hamiltonian_;
+    }
 
     void run() override;
 

tests/CMakeLists.txt

@@ -246,6 +246,8 @@ add_executable(testDensityMatrix
                ${CMAKE_SOURCE_DIR}/src/tools/mgmol_mpi_tools.cc
                ${CMAKE_SOURCE_DIR}/src/tools/random.cc
                ${CMAKE_SOURCE_DIR}/tests/ut_magma_main.cc)
+add_executable(testEnergyAndForces
+               ${CMAKE_SOURCE_DIR}/tests/EnergyAndForces/testEnergyAndForces.cc)
 
 if(${MAGMA_FOUND})
   add_executable(testOpenmpOffload
@@ -334,6 +336,14 @@ add_test(NAME testGramMatrix
 add_test(NAME testDensityMatrix
          COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS}
                  ${CMAKE_CURRENT_BINARY_DIR}/testDensityMatrix)
+add_test(NAME testEnergyAndForces
+         COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/EnergyAndForces/test.py
+         ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS}
+         ${CMAKE_CURRENT_BINARY_DIR}/testEnergyAndForces
+         ${CMAKE_CURRENT_SOURCE_DIR}/EnergyAndForces/mgmol.cfg
+         ${CMAKE_CURRENT_SOURCE_DIR}/EnergyAndForces/coords.in
+         ${CMAKE_CURRENT_SOURCE_DIR}/EnergyAndForces/lrs.in
+         ${CMAKE_CURRENT_SOURCE_DIR}/../potentials)
 
 if(${MAGMA_FOUND}) 
   add_test(NAME testOpenmpOffload
@@ -518,6 +528,7 @@ target_link_libraries(testBlacsContext PRIVATE ${SCALAPACK_LIBRARIES}
   ${BLAS_LIBRARIES} MPI::MPI_CXX)
 target_link_libraries(testSuperSampling PRIVATE MPI::MPI_CXX)
 target_link_libraries(testDirectionalReduce PRIVATE MPI::MPI_CXX)
+target_link_libraries(testEnergyAndForces PRIVATE mgmol_src)
 
 if(${MAGMA_FOUND})
   target_link_libraries(testDistVector PRIVATE ${SCALAPACK_LIBRARIES}

tests/EnergyAndForces/coords.in

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

tests/EnergyAndForces/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/EnergyAndForces/mgmol.cfg

@@ -0,0 +1,33 @@
+verbosity=1
+xcFunctional=LDA
+FDtype=Mehrstellen
+[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-9
+step_length=2.
+ortho_freq=10
+[Orbitals]
+initial_type=Gaussian
+initial_width=1.5
+temperature=10.
+nempty=1
+[Restart]
+output_level=0
+[DensityMatrix]
+mixing=0.5

tests/EnergyAndForces/test.py

@@ -0,0 +1,87 @@
+#!/usr/bin/env python
+import sys
+import os
+import subprocess
+import string
+
+print("Test EnergyAndForces...")
+
+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)
+
+flag=0
+forces=[]
+for line in lines:
+  if flag>0:
+    print(line)
+    words=line.split(b'    ')
+    forces.append(words[1].decode())
+    forces.append(words[2].decode())
+    forces.append(words[3].decode())
+    flag=flag-1
+  if line.count(b'Forces:'):
+    flag=2
+
+
+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("Check forces...")
+print(forces)
+flag=0
+for i in range(6):
+  diff=eval(forces[i+6])-eval(forces[i])
+  print(diff)
+  if abs(diff)>1.e-3:
+    print("Forces difference larger than tol")
+    flag=1
+if flag>0:
+  sys.exit(1)
+
+print("Test SUCCESSFUL!")
+sys.exit(0)

tests/EnergyAndForces/testEnergyAndForces.cc

@@ -0,0 +1,198 @@
+// Copyright (c) 2017, Lawrence Livermore National Security, LLC and
+// UT-Battelle, LLC.
+// Produced at the Lawrence Livermore National Laboratory and the Oak Ridge
+// National Laboratory.
+// LLNL-CODE-743438
+// All rights reserved.
+// This file is part of MGmol. For details, see https://github.com/llnl/mgmol.
+// Please also read this link https://github.com/llnl/mgmol/LICENSE
+
+#include "Control.h"
+#include "ExtendedGridOrbitals.h"
+#include "LocGridOrbitals.h"
+#include "MGmol.h"
+#include "MGmol_MPI.h"
+#include "MPIdata.h"
+#include "mgmol_run.h"
+
+#include <cassert>
+#include <iostream>
+#include <time.h>
+#include <vector>
+
+#include <boost/program_options.hpp>
+namespace po = boost::program_options;
+
+int main(int argc, char** argv)
+{
+    int mpirc = MPI_Init(&argc, &argv);
+    if (mpirc != MPI_SUCCESS)
+    {
+        std::cerr << "MPI Initialization failed!!!" << std::endl;
+        MPI_Abort(MPI_COMM_WORLD, 0);
+    }
+
+    MPI_Comm comm = MPI_COMM_WORLD;
+
+    /*
+     * Initialize general things, like magma, openmp, IO, ...
+     */
+    mgmol_init(comm);
+
+    /*
+     * read runtime parameters
+     */
+    std::string input_filename("");
+    std::string lrs_filename;
+    std::string constraints_filename("");
+
+    float total_spin = 0.;
+    bool with_spin   = false;
+
+    po::variables_map vm;
+
+    // read from PE0 only
+    if (MPIdata::onpe0)
+    {
+        read_config(argc, argv, vm, input_filename, lrs_filename,
+            constraints_filename, total_spin, with_spin);
+    }
+
+    MGmol_MPI::setup(comm, std::cout, with_spin);
+    MGmol_MPI& mmpi      = *(MGmol_MPI::instance());
+    MPI_Comm global_comm = mmpi.commGlobal();
+
+    /*
+     * Setup control struct with run time parameters
+     */
+    Control::setup(global_comm, with_spin, total_spin);
+    Control& ct = *(Control::instance());
+
+    ct.setOptions(vm);
+
+    int ret = ct.checkOptions();
+    if (ret < 0) return ret;
+
+    mmpi.bcastGlobal(input_filename);
+    mmpi.bcastGlobal(lrs_filename);
+
+    // Enter main scope
+    {
+        if (MPIdata::onpe0)
+        {
+            std::cout << "-------------------------" << std::endl;
+            std::cout << "Construct MGmol object..." << std::endl;
+            std::cout << "-------------------------" << std::endl;
+        }
+
+        MGmolInterface* mgmol;
+        if (ct.isLocMode())
+            mgmol = new MGmol<LocGridOrbitals>(global_comm, *MPIdata::sout,
+                input_filename, lrs_filename, constraints_filename);
+        else
+            mgmol = new MGmol<ExtendedGridOrbitals>(global_comm, *MPIdata::sout,
+                input_filename, lrs_filename, constraints_filename);
+
+        if (MPIdata::onpe0)
+        {
+            std::cout << "-------------------------" << std::endl;
+            std::cout << "MGmol setup..." << std::endl;
+            std::cout << "-------------------------" << std::endl;
+        }
+        mgmol->setup();
+
+        if (MPIdata::onpe0)
+        {
+            std::cout << "-------------------------" << std::endl;
+            std::cout << "Setup done..." << std::endl;
+            std::cout << "-------------------------" << std::endl;
+        }
+
+        // here we just use the atomic positions read in and used
+        // to initialize MGmol
+        std::vector<double> positions;
+        mgmol->getAtomicPositions(positions);
+        std::vector<short> anumbers;
+        mgmol->getAtomicNumbers(anumbers);
+        if (MPIdata::onpe0)
+        {
+            std::cout << "Positions:" << std::endl;
+            std::vector<short>::iterator ita = anumbers.begin();
+            for (std::vector<double>::iterator it = positions.begin();
+                 it != positions.end(); it += 3)
+            {
+                std::cout << *ita;
+                for (int i = 0; i < 3; i++)
+                    std::cout << "    " << *(it + i);
+                std::cout << std::endl;
+                ita++;
+            }
+        }
+
+        // compute energy and forces using all MPI tasks
+        // expect positions to be replicated on all MPI tasks
+        std::vector<double> forces;
+        double eks
+            = mgmol->evaluateEnergyAndForces(positions, anumbers, forces);
+
+        // print out results
+        if (MPIdata::onpe0)
+        {
+            std::cout << "Eks: " << eks << std::endl;
+            std::cout << "Forces:" << std::endl;
+            for (std::vector<double>::iterator it = forces.begin();
+                 it != forces.end(); it += 3)
+            {
+                for (int i = 0; i < 3; i++)
+                    std::cout << "    " << *(it + i);
+                std::cout << std::endl;
+            }
+        }
+
+        // move atoms one mesh spacing in all directions
+        Mesh* mymesh           = Mesh::instance();
+        const pb::Grid& mygrid = mymesh->grid();
+        const double hspacing[3]
+            = { mygrid.hgrid(0), mygrid.hgrid(1), mygrid.hgrid(2) };
+        int i = 0;
+        for (auto& pos : positions)
+        {
+            pos += hspacing[i % 3];
+            i++;
+        }
+
+        // evaluate energy and forces again
+        eks = mgmol->evaluateEnergyAndForces(positions, anumbers, forces);
+
+        // print out results
+        if (MPIdata::onpe0)
+        {
+            std::cout << "Eks: " << eks << std::endl;
+            std::cout << "Forces:" << std::endl;
+            for (std::vector<double>::iterator it = forces.begin();
+                 it != forces.end(); it += 3)
+            {
+                for (int i = 0; i < 3; i++)
+                    std::cout << "    " << *(it + i);
+                std::cout << std::endl;
+            }
+        }
+
+        delete mgmol;
+
+    } // close main scope
+
+    mgmol_finalize();
+
+    mpirc = MPI_Finalize();
+    if (mpirc != MPI_SUCCESS)
+    {
+        std::cerr << "MPI Finalize failed!!!" << std::endl;
+    }
+
+    time_t tt;
+    time(&tt);
+    if (onpe0) std::cout << " Run ended at " << ctime(&tt) << std::endl;
+
+    return 0;
+}