Training on stresses

[DanielMarchand]

Is your feature request related to a problem? Please describe.
DFT calculations usually include the stress output which can help out in training (if the settings are sufficiently converged). I would really like to use this for my calculations. I'm willing to invest of a bit of time implementing this but want to know if it's feasible.

[singraber]

Hello Daniel!

This is certainly not a trivial task but overall I would say that it is feasible. The expressions for the stress tensor can be found in this publication. An implementation would require two major parts:

1. Prediction As far as I can see it should be straightforward to implement the prediction for an existing NNP. All quantities necessary (forces, distances) are already computed during the force calculation and can afterwards be retrieved to compute also the stress tensor.

2. Training If I didn't miss anything the virial stress expression in the above paper would mean for the training that the derivatives of stress with respect to NN weights are just constructed from the force derivatives, hence this should not be very complicated either.

It's good that recently the Training class was further modularized to simplify the addition of new training "properties", such as stresses. This was done with the 4G-HDNNPs in mind (#77) where charges need to be fitted, but it will certainly be helpful here as well (most of these changes are already in the master but we may need to cherry-pick a few things from the 4G-HDNNP branch).

So yes, I think it can be done and if you have some experience in C/C++ coding it will probably not be very difficult. I will try to be helpful and point you to the corresponding parts in the code. However, I have to mention that the amount of time that I can spend on n2p2 development at this point is rather limited, considering that there are already other PRs in progress (#77, #92). So the major part of coding will have to rest on your shoulders.

Still, it would be great if you want to try it.. I would suggest to start a new feature branch and immediately create a "draft" pull request, so we can discuss the further proceeding there.

Best,

Andreas

[DanielMarchand]

Hi Andreas,

Thanks for the response!

I'm definitely interested. I may have to drop it if it becomes too difficult as it's not a key priority for my project.

For now, let's focus on (1). This will help me get a bit more familiar with the code and make sure I have the notation/units worked out correctly.

Let me dig into this a bit more and I'll message you when I'm ready. (I'll need to set up my development environment a bit more).

Daniel

[singraber]

Great, I am looking forward to this! Don't worry.. if you have to drop it I will still be useful as I will most likely continue it at some point.. it's definitely on my list but also not at a high priority.

Just to get you started here is some basic thing which should be added:

So far the configuration file format does not support stresses to be read in. Similarly, the Structure class has no member to store stresses either. Hence, there should be a new member array double stress[6] (and for training also double stressRef[6]) somewhere around here:

n2p2/src/libnnp/Structure.h

Lines 75 to 84 in 3f03449

	/// Potential energy determined by neural network.
	double energy;
	/// Reference potential energy.
	double energyRef;
	/// Charge determined by neural network potential.
	double charge;
	/// Reference charge.
	double chargeRef;
	/// Simulation box volume.
	double volume;

To support reading in the reference stresses I suggest to introduce a new keyword stress in the configuration file format, so users can add them in this way:

begin
...
stress <s_xx> <s_yy> <s_zz> <s_xy> <s_xz> <s_yz>
...
end

This can be simply done by adding another

else if (splitLine.at(0) == "stress")
{
    stressRef[0] = ...
    ...
}

section here:

n2p2/src/libnnp/Structure.cpp

Lines 131 to 243 in 3f03449

	void Structure::readFromLines(vector<string> const& lines)
	{
	size_t iBoxVector = 0;
	vector<string> splitLine;
	// read first line, should be keyword "begin".
	splitLine = split(reduce(lines.at(0)));
	if (splitLine.at(0) != "begin")
	{
	throw runtime_error("ERROR: Unexpected line content, expected"
	" \"begin\" keyword.\n");
	}
	for (vector<string>::const_iterator line = lines.begin();
	line != lines.end(); ++line)
	{
	splitLine = split(reduce(*line));
	if (splitLine.at(0) == "begin")
	{
	}
	else if (splitLine.at(0) == "comment")
	{
	size_t position = line->find("comment");
	string tmpLine = *line;
	comment = tmpLine.erase(position, splitLine.at(0).length() + 1);
	}
	else if (splitLine.at(0) == "lattice")
	{
	if (iBoxVector > 2)
	{
	throw runtime_error("ERROR: Too many box vectors.\n");
	}
	box[iBoxVector][0] = atof(splitLine.at(1).c_str());
	box[iBoxVector][1] = atof(splitLine.at(2).c_str());
	box[iBoxVector][2] = atof(splitLine.at(3).c_str());
	iBoxVector++;
	if (iBoxVector == 3)
	{
	isPeriodic = true;
	if (box[0][1] > numeric_limits<double>::min() ||
	box[0][2] > numeric_limits<double>::min() ||
	box[1][0] > numeric_limits<double>::min() ||
	box[1][2] > numeric_limits<double>::min() ||
	box[2][0] > numeric_limits<double>::min() ||
	box[2][1] > numeric_limits<double>::min())
	{
	isTriclinic = true;
	}
	calculateInverseBox();
	calculateVolume();
	}
	}
	else if (splitLine.at(0) == "atom")
	{
	atoms.push_back(Atom());
	atoms.back().index = numAtoms;
	atoms.back().indexStructure = index;
	atoms.back().tag = numAtoms; // Implicit conversion!
	atoms.back().r[0] = atof(splitLine.at(1).c_str());
	atoms.back().r[1] = atof(splitLine.at(2).c_str());
	atoms.back().r[2] = atof(splitLine.at(3).c_str());
	atoms.back().element = elementMap[splitLine.at(4)];
	atoms.back().chargeRef = atof(splitLine.at(5).c_str());
	atoms.back().fRef[0] = atof(splitLine.at(7).c_str());
	atoms.back().fRef[1] = atof(splitLine.at(8).c_str());
	atoms.back().fRef[2] = atof(splitLine.at(9).c_str());
	atoms.back().numNeighborsPerElement.resize(numElements, 0);
	numAtoms++;
	numAtomsPerElement[elementMap[splitLine.at(4)]]++;
	}
	else if (splitLine.at(0) == "energy")
	{
	energyRef = atof(splitLine[1].c_str());
	}
	else if (splitLine.at(0) == "charge")
	{
	chargeRef = atof(splitLine[1].c_str());
	}
	else if (splitLine.at(0) == "end")
	{
	if (!(iBoxVector == 0 || iBoxVector == 3))
	{
	throw runtime_error("ERROR: Strange number of box vectors.\n");
	}
	break;
	}
	else
	{
	throw runtime_error("ERROR: Unexpected file content, "
	"unknown keyword \"" + splitLine.at(0) +
	"\".\n");
	}
	}
	for (size_t i = 0; i < numElements; i++)
	{
	if (numAtomsPerElement[i] > 0)
	{
	numElementsPresent++;
	}
	}
	if (isPeriodic)
	{
	for (vector<Atom>::iterator it = atoms.begin();
	it != atoms.end(); ++it)
	{
	remap((*it));
	}
	}
	return;
	}

Best,

Andi

[DanielMarchand]

Hi Andi,

Please find a pull request #103. I've added some notes for discussion.

By the way I can't seem to build the code on anything other than the main branch. It seems like sed complains about an empty variable. Can you confirm it's possible to build the code outside of master?

Best,

Daniel

[singraber]

Hi Daniel,

I'm afraid the 4G-HDNNP-prediction branch is currently broken (which can be seen in the CI log). It's under heavy development and I am aware that there are a few things to fix...

However, an issue with sed is new to me, can you post the error maybe? Using sed is always a bit problematic in makefiles since there are different feature sets on OSes.. are you working on a Mac?

Best,

Andi

[DanielMarchand]

Hi Andi,

Actually, I had made a branch from master (perhaps I should change this to the 4G-HDNNP-prediction branch). Turned out it was my branch naming convention that was causing the crash. I've made a PR with a fix #106. I often use forward slashes in my branch name: feat/my_feature, patch/my_patch, ...

My PR will still break if someone pushes a branch with a plus sign (+) but I think that should be comparatively rare. The ideal fix would be simply to strip all special symbols first.