The unstable Metadynamics simulation using lammps

[wxsongsh]

Dear SSAGES Developers,

I added a new CV (Steinhardt order parameter) to SSAGES and carried out a Metadynamics simulation using lammps. Using NPT ensemble, the simulation is not stable and will break after many steps (less than 100000 steps). The system temperature always increases before the simulation break.

The output CV value is right by comparing the calculated value with the standard values of many crystalline phases. But I am not sure whether the deduced force is correct. The following is the Error imformation:

Step TotEng PotEng KinEng Press **Temp** Volume Lx Ly Lz
       0   -3366.7468   -3481.8529    115.10616    7472.1931         **1300**    11872.637       22.813       22.813       22.813
    1000    -3242.864   -3361.2629    118.39894   -4372.1207    **1337.1884**    12357.403    23.119356    23.119356    23.119356
...

  476000   -3144.3593   -3259.9841    115.62484   -1306.9965     **1305.858**    12614.993    23.278894    23.278894    23.278894
  477000    -3133.837   -3253.7667    119.92967    5327.3924    **1354.4764**     12575.77    23.254742    23.254742    23.254742
  478000   -3129.5073   -3250.0644    120.55704    4111.6205    **1361.5618**     12584.54    23.260147    23.260147    23.260147
  479000   -3145.7521    -3261.559    115.80683   -1659.5989    **1307.9133**    12623.331    23.284022    23.284022    23.284022
  480000   -3095.6606   -3221.6845    126.02384     4605.021    **1423.3034**    12776.289    23.377689    23.377689    23.377689
  481000   -2966.9635   -3146.1461    179.18263   -12047.642    **2023.6747**    13529.082    23.828102    23.828102    23.828102
  482000   -2991.1672   -3153.9548    162.78762    -9805.879    **1838.5107**    13373.685    23.736519    23.736519    23.736519
  483000    -2957.497   -3143.1443    185.64723    11727.174    **2096.6853**    13106.295    23.577259    23.577259    23.577259
ERROR: Non-numeric pressure - simulation unstable (../fix_nh.cpp:1038)
Last command:  run             1000000

Contents of Meta.json:

{
    "walkers" : 1,
    "input" : "in.metal",
    "CVs": [
        {
            "qatom": false,
            "qvmean": true,
            "short_cutoff": 2.8,
            "long_cutoff": 3.6,
            "q_dimension": 6,
            "type": "Q_Steinhardt"
        }
     ],
    "methods" : [
        {
            "type" : "Metadynamics",
            "widths" : [
                0.02
            ],
            "height" : 0.556075, (units in eV)
            "lower_bounds" : [-0.1],
            "upper_bounds" : [1.1],
            "lower_bound_restraints" : [0],
            "upper_bound_restraints" : [0],
            "hill_frequency" : 500
        }
    ]
}

In the CVXX.cpp file, I provide the value and grad vector values for each processor. What's wrong in my code? Is it the wrong bias force added?

Thanks! Looking forwards to your reply.

[hsidky]

If your CV is stable in the NVT ensemble, but not the NPT ensemble, then it's likely due to the fact that you have not implemented the CV "box gradient." This is the contribution of the CV to the virial. It is necessary for NPT ensembles.

You will need to calculate and set boxgrad_ in your CV. You can see an example of this for the ParticleSeparationCV here: https://github.com/SSAGESproject/SSAGES/blob/b0e010ce2115473bb72c681cd85abd4b5f33ca43/src/CVs/ParticleSeparationCV.h#L159

Hope this helps!

[mquevill]

The boxgrad_ can be calculated with either relative positions (rij in Particle SeparationCV is the vector between two points) or from absolute positions (using pos[i], shown below).

Label idx;
snapshot.GetLocalIndices(atomids_, &idx);
boxgrad_ = Matrix3::Zero();
for(auto& i : idx)
{
    boxgrad_ += grad_[i]*pos[i].transpose();
}

We have had discussions about standardizing this calculation, since the above code is generalizable to most CVs.

[wxsongsh]

Thanks very much! After adding the boxgrad_, the simulation is stable now. I do it like this:

// Initialize gradient.
std::fill(grad_.begin(), grad_.end(), Vector3{0,0,0});
grad_.resize(n, Vector3{0,0,0});
boxgrad_ = Matrix3::Zero();

for(int i=0;i<n;i++)
{
    int id=atomid[i]-1;
    grad_[i][0]=bias[id][0];//bias_force is a na(system total atom number)*3 array.
    grad_[i][1]=bias[id][1];
    grad_[i][2]=bias[id][2];

    boxgrad_ += grad_[i]*pos[i].transpose();
}

I am confused why boxgrad as well as grad is not gathered using the following sentence:

MPI_Allreduce(MPI_IN_PLACE,&boxgrad_,3,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);

[mquevill]

I'm glad to hear that your simulations are more stable with the added changes!

When a simulation is run with MPI, the atoms may be scattered across several processors. For calculations that require values from all processors, we have to do a collective operation, like MPI_Allreduce. However, the gradient in grad_ will only be affecting the atom for which it is calculated. Therefore, we do not have to reduce those values.

In your snippet of code, you are using atomid[i]. This is not recommended, since it may break when running with more than 1 MPI process. We recommend using the snapshot.GetLocalIndices() function to index matrices and vectors. If you have 4 atoms on 2 MPI processes, it could be split up in many ways, and it could even be an uneven split. If you are only running with 1 MPI process, your method of indexing may work fine.

I'm not sure exactly how your bias is calculated, so you would need to adapt your code to this example. We recommend the following:

Label idx; // Label is std::vector<unsigned int>
snapshot.GetLocalIndices(atomids_, &idx); // Finds local indices on processor

n = snapshot.GetNumAtoms();
bias_force = std::vector<Vector3>;
std::fill(bias_force.begin(), bias_force.end(), Vector3{0,0,0});
bias_force.resize(n, Vector3{0,0,0});

// Loop over local processor indices and calculate bias
for(auto& i : idx)
{
    bias_force[i] = ... // However you calculate this
}

// If you need to do operations with reduced data
// Not necessary if all bias_force[i] are separate
MPI_Allreduce(MPI_IN_PLACE, bias_force.data(), bias_force.size(), MPI_DOUBLE, MPI_SUM, snapshot.GetCommunicator());

// Loop over local processor to save gradient
for(auto& i : idx)
{
    grad_[i] = bias_force[i];
    boxgrad_ += grad_[i]*pos[i].transpose();
}

[wxsongsh]

Thanks, mquevill. I see that the grad and boxgrad variables should not be reduced in the CVXX.cpp, which will be done by SSAGES in other place.

As to the atomid[i] in my added CV, each atom (as a center atom) and its neighbor atoms are considered together to analyze the local structure characteristics. I need add the bias force in the the neighbor atoms, and then the sum of the minus force for the center atom. In this scenario, the neighbor atoms are probably in other processors, so the bias_force variable should be a global one.

At first, I obtain the positions of all atoms (na*3 matrix) by adding a new function GetAllPositions in the Snapshot.cpp and fix_ssages.cpp files. Then, in each processor, I can get the neighbor atoms of each atom in the processor. Next, I analyze the bias force for the center atoms and its neighbor atoms. So, the bias_force variable should be defined as a na*3 matrix here. Finally, I reduce the na*3 bias_force like this:

MPI_Allreduce(MPI_IN_PLACE,&bias_force[0][0],3*na,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);//in my definition, bias_force is a na*3 arrary, not a vector.

In oder to assign the values of grad_, I do it in the following:

auto atomid=snapshot.GetAtomIDs();

for(int i=0;i<n;i++)
{
    int id=atomid[i]-1;
    grad_[i][0]=bias_force[id][0];//bias_force is a na*3 array.
    grad_[i][1]=bias_force[id][1];
    grad_[i][2]=bias_force[id][2];
}