testLJ.cpp

An example demonstrating the computation of Piola and Cauchy stress tensors using hardy and virial spherical averaging domains, and 1D and 2D grids

1. Open the configuration file (*.data in MDStressLab format) and read the number of particles

       int numberOfParticles;
       int referenceAndFinal= true; // Does the input file include reference and final configurations
       std::string configFileName= "config.data";
       std::ifstream file(configFileName);
       if(!file) MY_ERROR("ERROR: config.data could not be opened for reading!\n");
       file >> numberOfParticles;
       if (numberOfParticles < 0) MY_ERROR("Error: Negative number of particles.\n");

2. Initialize the BoxConfiguration and read the reference and current atomic coordinates from the configuration file. The reference configuration is an fcc Ar crystal in the relaxed state. The deformed/current configuration is the reference crystal strained in the \(y\)-direction.

       BoxConfiguration body{numberOfParticles,referenceAndFinal};
       body.read(configFileName,referenceAndFinal);

3. Initialize the Kim model

       std::string modelname= "LJ_Smoothed_Bernardes_1958_Ar__MO_764178710049_001";
       Kim kim(modelname);

4. Create 1D and 2D grids

       int ngrid=200;
       Grid<Current> current2DGrid(Vector3d(0,0,0),Vector3d(60,60,60),ngrid,ngrid);
       Grid<Reference> reference2DGrid(Vector3d(0,0,0),Vector3d(60,60,60),ngrid,ngrid);
    
       ngrid= 500;
       Grid<Reference> reference1DGrid(Vector3d(0,0,0),Vector3d(60,60,60),ngrid);

5. Construct hardy and virial MethodSphere objects

       MethodSphere hardy5(5.0,"hardy");
       MethodSphere virial5(5.0,"virial");

6. Construct Stress objects using the MethodSphere and Grid objects

       Stress<MethodSphere,Cauchy> hardyCauchy("hardyCauchy",hardy5,&current2DGrid);
       Stress<MethodSphere,Cauchy> virialCauchy("virialCauchy",virial5,&current2DGrid);
    
       Stress<MethodSphere,Piola> hardyPiola("hardyPiola",hardy5,&reference1DGrid);
       Stress<MethodSphere,Piola> virialPiola("virialPiola",virial5,&reference1DGrid);

7. Calculate stresses. The following snippet shows that stresses can be calculated all at once or with various combinations.

       calculateStress(body,kim,
                       std::tie(),
                       std::tie());
    
       calculateStress(body,kim,
                       std::tie(hardyCauchy));
    
       calculateStress(body,kim,
                       std::tie(virialPiola));
    
       calculateStress(body,kim,
                       std::tie(hardyPiola,virialPiola),
                       std::tie(hardyCauchy,virialCauchy));
    

8. Write stresses

       hardyCauchy.write();
       hardyPiola.write();
       virialCauchy.write();
       virialPiola.write();

9. We compare our results with the exact results for unit testing purposes.

       compareStress("hardyCauchy");
       compareStress("hardyPiola");
       compareStress("virialCauchy");
       compareStress("virialPiola");

Visual comparison

• See Visualization Utilities for contour plotting

  Hardy Cauchy stress [View high‑res PDF]

  Virial Cauchy stress [View high‑res PDF]

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 "../compareStress.cpp"
 
int main()
{
    int numberOfParticles;
    int referenceAndFinal= true; // Does the input file include reference and final configurations
    std::string configFileName= "config.data";
    std::ifstream file(configFileName);
    if(!file) MY_ERROR("ERROR: config.data could not be opened for reading!\n");
    file >> numberOfParticles;
    if (numberOfParticles < 0) MY_ERROR("Error: Negative number of particles.\n");
    BoxConfiguration body{numberOfParticles,referenceAndFinal};
    body.read(configFileName,referenceAndFinal);
    std::string modelname= "LJ_Smoothed_Bernardes_1958_Ar__MO_764178710049_001";
    Kim kim(modelname);
    int ngrid=200;
    Grid<Current> current2DGrid(Vector3d(0,0,0),Vector3d(60,60,60),ngrid,ngrid);
    Grid<Reference> reference2DGrid(Vector3d(0,0,0),Vector3d(60,60,60),ngrid,ngrid);
 
    ngrid= 500;
    Grid<Reference> reference1DGrid(Vector3d(0,0,0),Vector3d(60,60,60),ngrid);
    MethodSphere hardy5(5.0,"hardy");
    MethodSphere virial5(5.0,"virial");
    Stress<MethodSphere,Cauchy> hardyCauchy("hardyCauchy",hardy5,&current2DGrid);
    Stress<MethodSphere,Cauchy> virialCauchy("virialCauchy",virial5,&current2DGrid);
 
    Stress<MethodSphere,Piola> hardyPiola("hardyPiola",hardy5,&reference1DGrid);
    Stress<MethodSphere,Piola> virialPiola("virialPiola",virial5,&reference1DGrid);
    calculateStress(body,kim,
                    std::tie(),
                    std::tie());
 
    calculateStress(body,kim,
                    std::tie(hardyCauchy));
 
    calculateStress(body,kim,
                    std::tie(virialPiola));
 
    calculateStress(body,kim,
                    std::tie(hardyPiola,virialPiola),
                    std::tie(hardyCauchy,virialCauchy));
 
    hardyCauchy.write();
    hardyPiola.write();
    virialCauchy.write();
    virialPiola.write();
    compareStress("hardyCauchy");
    compareStress("hardyPiola");
    compareStress("virialCauchy");
    compareStress("virialPiola");
    return 0;
}