crack/main.cpp

Compares Hardy stress computation using projected forces vs. process_dEdr with multiple OpenKIM models. The system is a single crystal Si with a crack under mode 1 loading

• Projected forces: Interatomic forces derived from forces.
• process_dEdr: Interatomic forces derived from analytical derivatives of energy w.r.t distances, provided by KIM models (if available).

For each model:

• We attempt to compute stress using projected forces first.
• Subsequently, we compute stress using process_dEdr if the model supports it.
• Both approaches use the MethodSphere class to compute Hardy stress on a spatial grid.

1. Query KIM models: For each model in the list, extract its lattice constant and elastic constants. In this demo, the list contains only one model.

       std::vector<std::string> modelnames={
               "Tersoff_LAMMPS_Tersoff_1988T3_Si__MO_186459956893_004"
               };
    
       std::cout << "Querying models." << std::endl;
       for (int i = 0; i < modelnames.size(); ++i) {
           std::cout << "Model: " << modelnames[i] << std::endl;
           std::regex rgx("MO_\\d+_\\d{3}");
    
           std::smatch match;
           if (!std::regex_search(modelnames[i], match, rgx))
               std::cout << "Regex error!" << std::endl;
    
    
           std::string latticeConstant= kim_query(
                   "get_lattice_constant_cubic",
                   {
                           "model=[\"" + match.str()+ "\"]",
                           "crystal=[\"diamond\"]",
                           "species=[\"Si\"]",
                           "units=[\"angstrom\"]"
                   }
           );
           std::cout << "Lattice constant = " << latticeConstant << std::endl;
           std::string elasticConstants= kim_query(
                   "get_elastic_constants_isothermal_cubic",
                   {
                           "model=[\"" + match.str()+ "\"]",
                           "crystal=[\"diamond\"]",
                           "species=[\"Si\"]",
                           "units=[\"eV/angstrom^3\"]"
                   }
           );
           std::cout << "Elastic constants = " << elasticConstants << std::endl;
           std::cout << "------------------------------------------ " << std::endl;
       }

2. Load atomic configuration: Read a MDStressLab-style configuration (config.data).

       int numberOfParticles;
       int referenceAndFinal= true;
       std::string configFileName= "config.data";
       std::ifstream file(configFileName);
       if(!file) MY_ERROR("ERROR: config.dat could not be opened for reading!");
    
       file >> numberOfParticles;
       if (numberOfParticles < 0) MY_ERROR("Error: Negative number of particles.\n");
    
       BoxConfiguration body{numberOfParticles,referenceAndFinal};
       body.read(configFileName,referenceAndFinal);

3. Load grid: Read grid points for stress evaluation from grid_cauchy.data.

       int ngrid;
       ngrid= 90601;
       Grid<Current> gridFromFile(ngrid);
       gridFromFile.read("grid_cauchy.data");

4. Compute and compare Hardy stress:
   • Start with stress computation using projected forces (always possible).
   • Try process_dEdr method, if available.
   • Write output stress fields for both methods (if available).

         MethodSphere hardy(6,"hardy");
         for (const auto modelname : modelnames)
         {
             Kim kim(modelname);
             try {
                 Stress<MethodSphere,Cauchy> hardyStress(hardy,&gridFromFile);
      
                 calculateStress(body, kim,
                                 std::tie(),
                                 std::tie(hardyStress), true);
                 hardyStress.write("project_hardy_" + modelname);
             }
             catch(const std::runtime_error& e){
                 std::cout << e.what() << std::endl;
                 std::cout << "Compute stress with projected forces failed. Moving on" << std::endl;
             }
      
             try{
                 Stress<MethodSphere,Cauchy> hardyStress(hardy,&gridFromFile);
      
                 // Calculate stress using the process_dedr, if possible
                 calculateStress(body, kim,
                                 std::tie(),
                                 std::tie(hardyStress));
                 hardyStress.write("hardy_" + modelname);
             }
             catch(const std::runtime_error& e){
                 std::cout << e.what() << std::endl;
                 std::cout << "Compute stress with process_dedr failed. Moving on" << std::endl;
             }
         }

Output:

For each KIM model, this example will generate:

• project_hardy_<modelname>: Stress computed using projected forces.
• hardy_<modelname>: Stress computed using process_dEdr.

Visual comparison

• See Visualization Utilities for contour plotting

  Projected force–based stress	process_dEdr–based stress
  Projected force–based stress	process_dEdr–based stress
  Projected force–based stress	process_dEdr–based stress

Full code:

/*
 * main.cpp
 *
 *  Created on: Nov 3, 2019
 *      Author: Nikhil Admal
 */
#include "MethodSphere.h"
#include <string>
#include <iostream>
#include <tuple>
#include <fstream>
#include "BoxConfiguration.h"
#include "calculateStress.h"
#include "Grid.h"
#include "typedef.h"
#include "kim_query.cpp"
#include <regex>
 
 
int main()
{
    std::vector<std::string> modelnames={
            "Tersoff_LAMMPS_Tersoff_1988T3_Si__MO_186459956893_004"
            };
 
    std::cout << "Querying models." << std::endl;
    for (int i = 0; i < modelnames.size(); ++i) {
        std::cout << "Model: " << modelnames[i] << std::endl;
        std::regex rgx("MO_\\d+_\\d{3}");
 
        std::smatch match;
        if (!std::regex_search(modelnames[i], match, rgx))
            std::cout << "Regex error!" << std::endl;
 
 
        std::string latticeConstant= kim_query(
                "get_lattice_constant_cubic",
                {
                        "model=[\"" + match.str()+ "\"]",
                        "crystal=[\"diamond\"]",
                        "species=[\"Si\"]",
                        "units=[\"angstrom\"]"
                }
        );
        std::cout << "Lattice constant = " << latticeConstant << std::endl;
        std::string elasticConstants= kim_query(
                "get_elastic_constants_isothermal_cubic",
                {
                        "model=[\"" + match.str()+ "\"]",
                        "crystal=[\"diamond\"]",
                        "species=[\"Si\"]",
                        "units=[\"eV/angstrom^3\"]"
                }
        );
        std::cout << "Elastic constants = " << elasticConstants << std::endl;
        std::cout << "------------------------------------------ " << std::endl;
    }
    int numberOfParticles;
    int referenceAndFinal= true;
    std::string configFileName= "config.data";
    std::ifstream file(configFileName);
    if(!file) MY_ERROR("ERROR: config.dat could not be opened for reading!");
 
    file >> numberOfParticles;
    if (numberOfParticles < 0) MY_ERROR("Error: Negative number of particles.\n");
 
    BoxConfiguration body{numberOfParticles,referenceAndFinal};
    body.read(configFileName,referenceAndFinal);
    int ngrid;
    ngrid= 90601;
    Grid<Current> gridFromFile(ngrid);
    gridFromFile.read("grid_cauchy.data");
    MethodSphere hardy(6,"hardy");
    for (const auto modelname : modelnames)
    {
        Kim kim(modelname);
        try {
            Stress<MethodSphere,Cauchy> hardyStress(hardy,&gridFromFile);
 
            calculateStress(body, kim,
                            std::tie(),
                            std::tie(hardyStress), true);
            hardyStress.write("project_hardy_" + modelname);
        }
        catch(const std::runtime_error& e){
            std::cout << e.what() << std::endl;
            std::cout << "Compute stress with projected forces failed. Moving on" << std::endl;
        }
 
        try{
            Stress<MethodSphere,Cauchy> hardyStress(hardy,&gridFromFile);
 
            // Calculate stress using the process_dedr, if possible
            calculateStress(body, kim,
                            std::tie(),
                            std::tie(hardyStress));
            hardyStress.write("hardy_" + modelname);
        }
        catch(const std::runtime_error& e){
            std::cout << e.what() << std::endl;
            std::cout << "Compute stress with process_dedr failed. Moving on" << std::endl;
        }
    }
    return 0;
}