Skip to content

Commit

Permalink
Charge neutrality check under PBC; cleanup
Browse files Browse the repository at this point in the history
  • Loading branch information
@philippmisof
philippmisof committed May 10, 2023
1 parent 15a2fd7 commit cf79881
Show file tree
Hide file tree
Showing 6 changed files with 44 additions and 74 deletions.
45 changes: 9 additions & 36 deletions src/libnnp/Mode.cpp
View file
Expand Up @@ -1634,8 +1634,7 @@ void Mode::calculateSymmetryFunctionGroups(Structure& structure,

void Mode::calculateAtomicNeuralNetworks(Structure& structure,
bool const derivatives,
string id,
bool const suppressOutput)
string id)
{
if (id == "") id = nnk.front();

Expand Down Expand Up @@ -1678,9 +1677,6 @@ void Mode::calculateAtomicNeuralNetworks(Structure& structure,
nn.getOutput(&(a.chi));
//double chi = a.chi;
if (normalize) a.chi = normalized("negativity", a.chi);
//if (!suppressOutput)
//log << strpr("Atom %5zu (%2s) chi: %16.8E\n",
// a.index, elementMap[a.element].c_str(), chi);
}
}
else if (id == "short")
Expand All @@ -1705,11 +1701,6 @@ void Mode::calculateAtomicNeuralNetworks(Structure& structure,
nn.calculateDEdG(&((a.dEdG).front()));
}
nn.getOutput(&(a.energy));
//if (!suppressOutput)
//log << strpr("Atom %5zu (%2s) energy: %16.8E\n",
// a.index + 1, elementMap[a.element].c_str(),
// a.energy
// + elements.at(a.element).getAtomicEnergyOffset());
}
}
}
Expand Down Expand Up @@ -1754,8 +1745,7 @@ void Mode::calculateAtomicNeuralNetworks(Structure& structure,

// TODO: Make this const?
void Mode::chargeEquilibration( Structure& structure,
bool const derivativesElec,
bool const suppressOutput)
bool const derivativesElec)
{
Structure& s = structure;

Expand All @@ -1782,17 +1772,13 @@ void Mode::chargeEquilibration( Structure& structure,
}
}

double const error = s.calculateElectrostaticEnergy(ewaldSetup,
hardness,
gammaSqrt2,
sigmaSqrtPi,
screeningFunction,
fourPiEps,
erfcBuf);


//if (!suppressOutput)
// log << strpr("Solve relative error: %16.8E\n", error);
s.calculateElectrostaticEnergy(ewaldSetup,
hardness,
gammaSqrt2,
sigmaSqrtPi,
screeningFunction,
fourPiEps,
erfcBuf);

// TODO: leave these 2 functions here or move it to e.g. forces? Needs to be
// executed after calculateElectrostaticEnergy.
Expand All @@ -1806,23 +1792,10 @@ void Mode::chargeEquilibration( Structure& structure,
fourPiEps);
}



//for (auto const& a : structure.atoms)
//{
// if (!suppressOutput)
// log << strpr("Atom %5zu (%2s) q: %16.8E\n",
// a.index, elementMap[a.element].c_str(), a.charge);
// //structure.charge += a.charge;
//}
//
//if (!suppressOutput)
//{
// log << strpr("Total charge: %16.8E (ref: %16.8E)\n",
// structure.charge, structure.chargeRef);

// log << strpr("Electrostatic energy: %16.8E\n", structure.energyElec);
//}
return;
}

Expand Down
11 changes: 2 additions & 9 deletions src/libnnp/Mode.h
View file
Expand Up @@ -336,28 +336,21 @@ class Mode
* (required for force calculation).
* @param[in] id Neural network ID to use. If empty, the first entry
* nnk.front() is used.
* @param[in] suppressOutput Turn on/off output generation (for example
* during training.
*/
// TODO: remove suppressOutput again
void calculateAtomicNeuralNetworks(
Structure& structure,
bool const derivatives,
std::string id = "",
bool const suppressOutput = false);
std::string id = "");
/** Perform global charge equilibration method.
*
* @param[in] structure Input structure.
* @param[in] derivativesElec Turn on/off calculation of dElecdQ and
* pElecpr (Typically needed for elecstrosttic
* forces).
* @param[in] suppressOutput Turn on/off output generation (for example
* during training.
*/
void chargeEquilibration(
Structure& structure,
bool const derivativesElec,
bool const suppressOutput = false);
bool const derivativesElec);
/** Calculate potential energy for a given structure.
*
* @param[in] structure Input structure.
Expand Down
9 changes: 9 additions & 0 deletions src/libnnp/Structure.cpp
View file
Expand Up @@ -238,6 +238,15 @@ void Structure::readFromLines(vector<string> const& lines)
{
throw runtime_error("ERROR: Strange number of box vectors.\n");
}
if (isPeriodic &&
abs(chargeRef) > 10*std::numeric_limits<double>::epsilon())
{
throw runtime_error("ERROR: In structure with index "
+ to_string(index)
+ "; if PBCs are applied, the\n"
"simulation cell has to be neutrally "
"charged.\n");
}
break;
}
else
Expand Down
8 changes: 3 additions & 5 deletions src/libnnpif/LAMMPS/InterfaceLammps.cpp
View file
Expand Up @@ -462,18 +462,16 @@ void InterfaceLammps::finalizeNeighborList()

void InterfaceLammps::process()
{
// TODO: remove this?
bool suppressOutput = true;
#ifdef N2P2_NO_SF_GROUPS
calculateSymmetryFunctions(structure, true);
#else
calculateSymmetryFunctionGroups(structure, true);
#endif
calculateAtomicNeuralNetworks(structure, true, "", suppressOutput);
calculateAtomicNeuralNetworks(structure, true, "");
if (nnpType == NNPType::HDNNP_4G)
{
chargeEquilibration(structure, true, suppressOutput);
calculateAtomicNeuralNetworks(structure, true, "short", suppressOutput);
chargeEquilibration(structure, true);
calculateAtomicNeuralNetworks(structure, true, "short");
ewaldSetup.logEwaldCutoffs(log, convLength * cflength);
}
calculateEnergy(structure);
Expand Down
39 changes: 19 additions & 20 deletions src/libnnptrain/Training.cpp
View file
Expand Up @@ -1429,25 +1429,25 @@ void Training::calculateError(
{
if ( stage == 1 )
{
calculateAtomicNeuralNetworks((*it), useForces, nnId, true);
chargeEquilibration((*it), false, true);
calculateAtomicNeuralNetworks((*it), useForces, nnId);
chargeEquilibration((*it), false);
}
else
{
if ( !it->hasCharges || (!it->hasAMatrix && useForces) )
{
calculateAtomicNeuralNetworks((*it), useForces,
"elec", true);
chargeEquilibration((*it), useForces, true);
"elec");
chargeEquilibration((*it), useForces);
}
calculateAtomicNeuralNetworks((*it), useForces, "short", true);
calculateAtomicNeuralNetworks((*it), useForces, "short");
calculateEnergy((*it));
if (useForces) calculateForces((*it));
}
}
else
{
calculateAtomicNeuralNetworks((*it), useForces, nnId, true);
calculateAtomicNeuralNetworks((*it), useForces, nnId);
calculateEnergy((*it));
if (useForces) calculateForces((*it));
}
Expand Down Expand Up @@ -2362,10 +2362,10 @@ void Training::update(string const& property)
{
if (!s.hasCharges)
{
calculateAtomicNeuralNetworks(s, derivatives, "elec", true);
chargeEquilibration(s, derivatives, true);
calculateAtomicNeuralNetworks(s, derivatives, "elec");
chargeEquilibration(s, derivatives);
}
calculateAtomicNeuralNetworks(s, derivatives, "short", true);
calculateAtomicNeuralNetworks(s, derivatives, "short");
calculateEnergy(s);
currentRmseFraction.at(b) = fabs(s.energyRef - s.energy)
/ (s.numAtoms
Expand Down Expand Up @@ -2395,10 +2395,10 @@ void Training::update(string const& property)
{
if (!s.hasAMatrix)
{
calculateAtomicNeuralNetworks(s, derivatives, "elec", true);
chargeEquilibration(s, derivatives, true);
calculateAtomicNeuralNetworks(s, derivatives, "elec");
chargeEquilibration(s, derivatives);
}
calculateAtomicNeuralNetworks(s, derivatives, "short", true);
calculateAtomicNeuralNetworks(s, derivatives, "short");
calculateForces(s);
Atom const& a = s.atoms.at(sC->a);
currentRmseFraction.at(b) =
Expand All @@ -2418,8 +2418,8 @@ void Training::update(string const& property)
// Assume stage 1.
if (nnpType == NNPType::HDNNP_4G)
{
calculateAtomicNeuralNetworks(s, derivatives,"",true);
chargeEquilibration(s, false, true);
calculateAtomicNeuralNetworks(s, derivatives, "");
chargeEquilibration(s, false);
Eigen::VectorXd QError;
double QErrorNorm;
calculateChargeErrorVec(s, QError, QErrorNorm);
Expand Down Expand Up @@ -2564,8 +2564,8 @@ void Training::update(string const& property)
{
if (nnpType == NNPType::HDNNP_4G && !s.hasCharges)
{
calculateAtomicNeuralNetworks(s, derivatives, "elec", true);
chargeEquilibration(s, derivatives, true);
calculateAtomicNeuralNetworks(s, derivatives, "elec");
chargeEquilibration(s, derivatives);
}
// Loop over atoms and calculate atomic energy contributions.
for (vector<Atom>::iterator it = s.atoms.begin();
Expand Down Expand Up @@ -2611,9 +2611,8 @@ void Training::update(string const& property)
{
if (!s.hasAMatrix)
{
calculateAtomicNeuralNetworks(s, derivatives, "elec",
true);
chargeEquilibration(s, derivatives, true);
calculateAtomicNeuralNetworks(s, derivatives, "elec");
chargeEquilibration(s, derivatives);
}
s.calculateDQdr(vector<size_t>{sC->a},
vector<size_t>{sC->c},
Expand Down Expand Up @@ -2718,7 +2717,7 @@ void Training::update(string const& property)
}

}
chargeEquilibration(s, false, true);
chargeEquilibration(s, false);

vector<Eigen::VectorXd> dQdChi;
s.calculateDQdChi(dQdChi);
Expand Down
6 changes: 2 additions & 4 deletions test/cpp/Example_nnp_train.h
View file
Expand Up @@ -41,7 +41,7 @@ void BoostDataContainer<Example_nnp_train>::setup()
e->rmseEnergyTest = 3.57910403E-04;
e->rmseForcesTrain = 2.46334310E-02;
e->rmseForcesTest = 1.22544569E-02;
e->accuracy = 2.5E-9;
e->accuracy = 1.0E-8;

examples.push_back(Example_nnp_train("H2O_RPBE-D3_norm-force"));
e = &(examples.back());
Expand All @@ -68,9 +68,7 @@ void BoostDataContainer<Example_nnp_train>::setup()
e->rmseEnergyTest = 1.30155783E-05;
e->rmseForcesTrain = 2.07286126E-04;
e->rmseForcesTest = 2.20371908E-04;
e->accuracy = 100.0 * std::numeric_limits<double>::epsilon();


e->accuracy = 1E-13;

return;
}
Expand Down

0 comments on commit cf79881

Please sign in to comment.