nnp-checkf.cpp File Reference

#include "Dataset.h"
#include "mpi-extra.h"
#include "utility.h"
#include <mpi.h>
#include <algorithm>
#include <cmath>
#include <cstdlib>
#include <fstream>
#include <iostream>
#include <map>
#include <string>
#include <vector>

Include dependency graph for nnp-checkf.cpp:

Go to the source code of this file.

Functions
int	main (int argc, char *argv[])

Function Documentation

main()

int main	(	int	argc,
		char *	argv[]
	)

Definition at line 33 of file nnp-checkf.cpp.

{
    int      numProcs = 0;
    int      myRank = 0;
    size_t   numOriginals = 0;
    double   delta = 0.0;
    ifstream dataFile;
    ofstream myLog;
    ofstream outFileForces;
    ofstream outFileSummary;
 
    if (argc > 2)
    {
        cout << "USAGE: " << argv[0] << " <<delta>>\n"
             << "       <<delta>> ... (optional) Displacement for central "
                "difference (default: 1.0e-4).\n"
             << "       Execute in directory with these NNP files present:\n"
             << "       - input.data (structure file)\n"
             << "       - input.nn (NNP settings)\n"
             << "       - scaling.data (symmetry function scaling data)\n"
             << "       - \"weights.%%03d.data\" (weights files)\n";
        return 1;
    }
 
    MPI_Init(&argc, &argv);
    MPI_Comm_size(MPI_COMM_WORLD, &numProcs);
    MPI_Comm_rank(MPI_COMM_WORLD, &myRank);
 
    if (argc == 2) delta = atof(argv[1]);
    else delta = 1.0E-4;
 
    Dataset dataset;
    if (myRank != 0) dataset.log.writeToStdout = false;
    myLog.open(strpr("nnp-checkf.log.%04d", myRank).c_str());
    dataset.log.registerStreamPointer(&myLog);
    dataset.setupMPI();
    dataset.initialize();
    dataset.loadSettingsFile();
    dataset.setupGeneric();
    bool normalize = dataset.useNormalization();
    dataset.setupSymmetryFunctionScaling();
    dataset.setupSymmetryFunctionStatistics(false, false, false, false);
    dataset.setupNeuralNetworkWeights();
 
    dataset.log << "\n";
    dataset.log << "*** ANALYTIC/NUMERIC FORCES CHECK *******"
                   "**************************************\n";
    dataset.log << "\n";
    dataset.log << strpr("Delta for symmetric difference quotient: %11.3E\n",
                         delta);
    if (myRank == 0)
    {
        string fileName = "checkf-forces.out";
        dataset.log << strpr("Individual analytic/numeric forces will be "
                             "written to \"%s\"\n", fileName.c_str());
        outFileForces.open(fileName.c_str());
        // File header.
        vector<string> title;
        vector<string> colName;
        vector<string> colInfo;
        vector<size_t> colSize;
        title.push_back(strpr("Comparison of analytic and numeric forces "
                              "(delta = %11.3E).", delta));
        colSize.push_back(10);
        colName.push_back("struct");
        colInfo.push_back("Structure index (starting with 1).");
        colSize.push_back(10);
        colName.push_back("atom");
        colInfo.push_back("Atom index (starting with 1).");
        colSize.push_back(3);
        colName.push_back("xyz");
        colInfo.push_back("Force component (0 = x, 1 = y, 2 = z).");
        colSize.push_back(24);
        colName.push_back("F_analytic");
        colInfo.push_back("Force computed from analytic derivative of NNP "
                          "energy.");
        colSize.push_back(24);
        colName.push_back("F_numeric");
        colInfo.push_back("Force computed numerically from symmetric "
                          "difference quotient.");
        appendLinesToFile(outFileForces,
                          createFileHeader(title, colSize, colName, colInfo));
 
        fileName = "checkf-summary.out";
        dataset.log << strpr("Per-structure summary of analytic/numeric force "
                             "comparison will be \n"
                             "written to \"%s\"\n",
                             fileName.c_str());
        outFileSummary.open(fileName.c_str());
        // File header.
        title.clear();
        colName.clear();
        colInfo.clear();
        colSize.clear();
        title.push_back(strpr("Per-structure summary of analytic vs. numeric "
                              "force comparison (delta = %11.3E).", delta));
        colSize.push_back(10);
        colName.push_back("struct");
        colInfo.push_back("Structure index (starting with 1).");
        colSize.push_back(16);
        colName.push_back("numForces");
        colInfo.push_back("Number of forces in this structure");
        colSize.push_back(16);
        colName.push_back("meanAbsError");
        colInfo.push_back("Mean over all absolute differences between "
                          "analytic and numeric forces in this structure.");
        colSize.push_back(16);
        colName.push_back("maxAbsError");
        colInfo.push_back("Maximum over all absolute differences between "
                          "analytic and numeric forces in this structure.");
        appendLinesToFile(outFileSummary,
                          createFileHeader(title, colSize, colName, colInfo));
    }
 
    // First, check how many originals are in "input.data".
    if (myRank == 0)
    {
        string fileName = "input.data";
        dataFile.open(fileName.c_str());
        numOriginals = dataset.getNumStructures(dataFile);
        dataset.log << strpr("Found %zu configurations in data file: %s.\n",
                             numOriginals,
                             fileName.c_str());
        dataFile.clear();
        dataFile.seekg(0);
    }
    MPI_Bcast(&numOriginals, 1, MPI_SIZE_T, 0, MPI_COMM_WORLD);
    dataset.log << strpr("Starting loop over %zu configurations...\n",
                          numOriginals);
 
    if (myRank == 0)
    {
        dataset.log << "\n";
        dataset.log << strpr("                      %10s %12s %12s  verdict\n",
                             "numForces", "meanAbsError", "maxAbsError");
        dataset.log << "-----------------------------------"
                       "--------------------------------\n";
    }
 
    bool warning = false;
    for (size_t is = 0; is < numOriginals; ++is)
    {
        // Read original, prepare and distribute modified copies.
        Structure original;
        original.setElementMap(dataset.elementMap);
        if (myRank == 0) original.readFromFile(dataFile);
        size_t numStructures = dataset.prepareNumericForces(original, delta);
        if (normalize) dataset.toNormalizedUnits();
 
        // Prepare arrays to collect central difference values and identifiers.
        vector<double> allEnergies(numStructures - 1, 0.0);
        vector<size_t> allAtoms(numStructures - 1, 0);
        vector<size_t> allXYZ(numStructures - 1, 0);
        vector<int> allSigns(numStructures - 1, 0);
 
        // Loop over all structures and compute symmetry functions and energies.
        for (auto& s : dataset.structures)
        {
            bool needForces = false;
            if (s.comment == "original") needForces = true;
            s.calculateNeighborList(dataset.getMaxCutoffRadius());
#ifdef N2P2_NO_SF_GROUPS
            dataset.calculateSymmetryFunctions(s, needForces);
#else
            dataset.calculateSymmetryFunctionGroups(s, needForces);
#endif
            dataset.calculateAtomicNeuralNetworks(s, needForces);
            dataset.calculateEnergy(s);
            if (needForces) dataset.calculateForces(s);
            s.freeAtoms(true);
        }
        if (normalize) dataset.toPhysicalUnits();
 
        // Loop over all structures collect energies and identifiers.
        for (auto& s : dataset.structures)
        {
            // Skip for original structure.
            if (s.comment == "original") continue;
            // Store central difference values in arrays (locally on each
            // processor, communicate later).
            vector<string> lsplit = split(s.comment);
            size_t count = atoi(lsplit.at(0).c_str());
            size_t iAtom = atoi(lsplit.at(1).c_str());
            size_t ixyz = atoi(lsplit.at(2).c_str());
            int sign = atoi(lsplit.at(3).c_str());
            allEnergies.at(count) = s.energy;
            allAtoms.at(count) = iAtom;
            allXYZ.at(count) = ixyz;
            allSigns.at(count) = sign;
        }
 
        // Collect data on rank 0.
        if (myRank == 0)
        {
            MPI_Reduce(MPI_IN_PLACE      , allEnergies.data(), numStructures - 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
            MPI_Reduce(MPI_IN_PLACE      , allAtoms.data()   , numStructures - 1, MPI_SIZE_T, MPI_SUM, 0, MPI_COMM_WORLD);
            MPI_Reduce(MPI_IN_PLACE      , allXYZ.data()     , numStructures - 1, MPI_SIZE_T, MPI_SUM, 0, MPI_COMM_WORLD);
            MPI_Reduce(MPI_IN_PLACE      , allSigns.data()   , numStructures - 1, MPI_INT   , MPI_SUM, 0, MPI_COMM_WORLD);
        }
        else
        {
            MPI_Reduce(allEnergies.data(), allEnergies.data(), numStructures - 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
            MPI_Reduce(allAtoms.data()   , allAtoms.data()   , numStructures - 1, MPI_SIZE_T, MPI_SUM, 0, MPI_COMM_WORLD);
            MPI_Reduce(allXYZ.data()     , allXYZ.data()     , numStructures - 1, MPI_SIZE_T, MPI_SUM, 0, MPI_COMM_WORLD);
            MPI_Reduce(allSigns.data()   , allSigns.data()   , numStructures - 1, MPI_INT   , MPI_SUM, 0, MPI_COMM_WORLD);
        }
 
        if (myRank == 0)
        {
            // Now temporarily store central difference energies in
            // fRef(+delta) and f(-delta) of original structure.
            for (size_t i = 0; i < numStructures - 1; ++i)
            {
                Atom& a = original.atoms.at(allAtoms.at(i));
                if (allSigns.at(i) == 1)
                {
                    a.fRef[allXYZ.at(i)] = allEnergies.at(i);
                }
                else if (allSigns.at(i) == -1)
                {
                    a.f[allXYZ.at(i)] = allEnergies.at(i);
                }
            }
 
            double maxAbsError = 0.0;
            double meanAbsError = 0.0;
            for (size_t i = 0; i < original.atoms.size(); ++i)
            {
                Atom const& ref = dataset.structures.at(0).atoms.at(i);
                Atom& a = original.atoms.at(i);
                for (size_t j = 0; j < 3; ++j)
                {
                    a.f[j] = - (a.fRef[j] - a.f[j]) / (2.0 * delta);
                    a.fRef[j] = ref.f[j];
                    double const error = fabs(a.f[j] - a.fRef[j]);
                    meanAbsError += error;
                    maxAbsError = max(error, maxAbsError);
                    outFileForces << strpr("%10zu %10zu %3zu %24.16E "
                                           "%24.16E\n",
                                           is + 1,
                                           i + 1,
                                           j,
                                           a.fRef[j],
                                           a.f[j]);
                }
            }
            size_t const numForces = 3 * original.atoms.size();
            meanAbsError /= numForces;
            dataset.log << strpr("Configuration %6zu: %10zu %12.3E %12.3E",
                                 is + 1, numForces, meanAbsError, maxAbsError);
            outFileSummary << strpr("%10zu %16zu %16.8E %16.8E\n",
                                    is + 1,
                                    numForces,
                                    meanAbsError,
                                    maxAbsError);
            if (maxAbsError > 10 * delta * delta)
            {
                dataset.log << "  WARNING!\n";
                warning = true;
            }
            else dataset.log << "  OK.\n";
        }
    }
 
    if (myRank == 0)
    {
        if (warning)
        {
            dataset.log << "\n"; 
            dataset.log << "IMPORTANT: Some warnings were issued. By default, this happens if the maximum\n"
                           "           absolute error (\"maxAbsError\") is higher than 10 * delta².\n"
                           "           However, this does NOT mean that analytic forces are incorrect!\n"
                           "           Repeat this analysis with a different delta and check whether\n"
                           "           the error scales with O(delta²). The prefactor for your system\n"
                           "           could be higher than 10, hence, as long as there is a O(delta²)\n"
                           "           scaling the analytic forces are probably correct.\n";
        }
        outFileForces.close();
        outFileSummary.close();
        dataFile.close();
    }
 
    dataset.log << "*****************************************"
                   "**************************************\n";
 
    myLog.close();
 
    MPI_Finalize();
 
    return 0;
}

References nnp::appendLinesToFile(), nnp::Structure::atoms, nnp::Mode::calculateAtomicNeuralNetworks(), nnp::Mode::calculateEnergy(), nnp::Mode::calculateForces(), nnp::Mode::calculateSymmetryFunctionGroups(), nnp::Mode::calculateSymmetryFunctions(), nnp::createFileHeader(), nnp::Mode::elementMap, nnp::Atom::f, nnp::Atom::fRef, nnp::Mode::getMaxCutoffRadius(), nnp::Dataset::getNumStructures(), nnp::Mode::initialize(), nnp::Mode::loadSettingsFile(), nnp::Mode::log, MPI_SIZE_T, nnp::Dataset::prepareNumericForces(), nnp::Structure::readFromFile(), nnp::Log::registerStreamPointer(), nnp::Structure::setElementMap(), nnp::Mode::setupGeneric(), nnp::Dataset::setupMPI(), nnp::Mode::setupNeuralNetworkWeights(), nnp::Mode::setupSymmetryFunctionScaling(), nnp::Mode::setupSymmetryFunctionStatistics(), nnp::split(), nnp::strpr(), nnp::Dataset::structures, nnp::Dataset::toNormalizedUnits(), nnp::Dataset::toPhysicalUnits(), nnp::Mode::useNormalization(), and nnp::Log::writeToStdout.

Here is the call graph for this function: