Pressure Minimisation Failing While Trying to optimise Lattice Parameters at 0 K

LAMMPS LAMMPS Beginners
1

I am working with the Rutile TiO2 structure, specifically a 4 x 4 x 2 bulk supercell, with periodic boundary conditions enforced in all three directions. My goal is to optimise the lattice parameters of this supercell at absolute zero temperature. Essentially, I want to relax the structure anisotropically to zero pressure (given Rutile’s tetragonal symmetry). I trained a DeepMD potential beforehand which I use to compute pairwise interactions.

I wrote a script reading in the data file with an initial set of atomic coordinates that I created from a CIF file obtained using the Materials Project Database. I then used fix box/relax with aniso 0.0 and vmax initially set to 1e-3, (I made vmax smaller (going to 1e-5) in a restarted run as well) then ran a minimisation of the simulation box followed by simultaneous refinement of atomic positions.

During the initial run, the pressure decreased monotonically for a large part of its run, then inexplicably flattened out and decreased at a much slower rate. When it terminated, the stopping criterion was “Energy Tolerance”, which makes me think that the algorithm swiftly converged to a local energy minimum and was unable to get out of it. In later runs, I would also get a “linesearch alpha is zero” stopping criterion.

Interestingly, I also tried this with a 5 x 5 x 5 supercell. There, I used the same approach with isotropic relaxation and it converged to zero pressure within a few hundred steps, although its local symmetry was significantly altered. Even isotropic relaxation, when applied to the 4 x 4 x 2 supercell, faces energy tolerance criterion stoppage issues at large non-zero pressures.

I also tried the approach taken by one of the official tutorials here. This approach, when applied to the 4 x 4 x 2 system, also stopped prematurely citing the same energy tolerance criterion, and I cannot figure out why this is happening.

I have included my input file, data file and log file from the initial run of the 4 x 4x 2 system (wherein I tried anisotropically relaxing the structure). I would appreciate any insights into why the pressure and forces are failing to converge. I am also happy to share additional scripts from other runs if need be.

Input Script:

########################################################################
units          metal
atom_style     atomic
dimension      3
boundary       p p p
read_data    rutile_442.lmp
#
pair_style     deepmd frozen_model.pb
pair_coeff     * *
neighbor       2.0 bin
neigh_modify   every 1 delay 0 check yes
#
thermo         1
thermo_style   custom step pe press lx ly lz
write_restart  restart.file
#
#########################################################################
## 1.  USER PARAMETERS FOR CONVERGENCE
#########################################################################
## Target and tolerance (bar)
variable Ptarget   equal 0.0
#
 # (A) Anisotropic box relaxation
# ######################################################################
  min_style cg
  fix     aniso all box/relax aniso ${Ptarget} vmax 1e-3
  minimize 1e-8 1e-8 10000 100000
#
#######################################################################
# (B) Final atomic positions refinement
 minimize 1e-8 1e-8 10000 100000

Data (containing atomic coordinates):

# Ti2 O4
  
         192  atoms
           2  atom types
 
      0.000000000000      18.399349280000  xlo xhi
      0.000000000000      18.399349280000  ylo yhi
      0.000000000000       5.918427120000  zlo zhi
 
Masses
 
            1   47.86700000             # Ti
            2   15.99900000             # O
 
Atoms # atomic
 
         1    1         2.299918660000       2.299918660000       0.000000000000
         2    1         0.000000000000       0.000000000000       1.479606780000
         3    2         3.200008800991       3.200008800991       1.479606780000
         4    2         0.900090140991       3.699747179009       0.000000000000
         5    2         3.699747179009       0.900090140991       0.000000000000
        ....
        (187 more rows)

Log file:

CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE

Generated 0 of 1 mixed pair_coeff terms from geometric mixing rule
Neighbor list info ...
  update: every = 1 steps, delay = 0 steps, check = yes
  max neighbors/atom: 2000, page size: 100000
  master list distance cutoff = 9
  ghost atom cutoff = 9
  binsize = 4.5, bins = 5 5 2
  1 neighbor lists, perpetual/occasional/extra = 1 0 0
  (1) pair deepmd, perpetual
      attributes: full, newton on
      pair build: full/bin/atomonly
      stencil: full/bin/3d
      bin: standard
WARNING: Energy due to 3 extra global DOFs will be included in minimizer energies
 (src/min.cpp:219)
Per MPI rank memory allocation (min/avg/max) = 4.301 | 4.301 | 4.301 Mbytes
   Step         PotEng         Press            Lx             Ly             Lz      
         0  -1515.9694      428834.72      18.399349      18.399349      5.9184271    
         1  -1517.3283      405322.28      18.417728      18.417749      5.9226245    
         2  -1518.6655      381601.34      18.435451      18.436148      5.9266257    
         3  -1519.9085      355670.67      18.453065      18.454547      5.9304581    
         4  -1521.0843      327694.44      18.470556      18.472947      5.9341003    
         5  -1522.1255      297607.04      18.487888      18.491346      5.9375196    
         6  -1523.0799      266619.99      18.505043      18.509745      5.9406847    
         7  -1523.8857      233321.48      18.521978      18.528145      5.9435459    
         8  -1524.5943      199807.89      18.538648      18.546544      5.9460411    
         9  -1525.2012      166346.18      18.554983      18.564943      5.9480827    
        10  -1525.7015      133455.7       18.570891      18.583343      5.9495536    
        11  -1526.1057      102181.25      18.586252      18.601742      5.9502957    
        12  -1526.4397      73456.605      18.600908      18.620141      5.9501022    
        13  -1526.7146      48425.305      18.614654      18.638541      5.9487085    
        14  -1526.9572      28360.847      18.627265      18.65694       5.9458194    
        15  -1527.1998      14276.187      18.638567      18.67534       5.9412092    
        16  -1527.4445      6596.3123      18.648         18.692658      5.9352908    
        17  -1527.6741      3525.0142      18.655239      18.70726       5.9293723    
        18  -1527.8903      1968.9798      18.661674      18.72088       5.9234539    
        19  -1528.1018      1170.3746      18.667711      18.734014      5.9175355    
        20  -1528.3106      726.97072      18.673541      18.746896      5.9116171    
        21  -1528.5176      463.72301      18.679252      18.759631      5.9056986    
        22  -1528.7232      293.15733      18.68489       18.772272      5.8997802    
        23  -1528.8341      97.661251      18.689945      18.783648      5.8944251    
        24  -1528.9794     -808.00065      18.708345      18.782988      5.8886721    
        25  -1529.0513     -1183.7272      18.725033      18.768131      5.8879541    
        26  -1530.8136      11939.991      18.743433      18.764985      5.8826662    
        27  -1531.1576      12352.859      18.758846      18.769248      5.8767478    
        28  -1531.3968      10544.744      18.769871      18.780264      5.8708294    
        29  -1531.4358      8087.7585      18.788105      18.777379      5.8670713    
        30  -1531.4389     -913.62874      18.79156       18.779916      5.8683218    
        31  -1531.4401     -875.05495      18.790627      18.779133      5.8688243    
        32  -1531.4414      534.77054      18.789038      18.777862      5.8692071    
        33  -1531.4425     -195.5915       18.788468      18.777536      5.869724     
        34  -1531.443      -105.45002      18.78797       18.777174      5.869952     
        35  -1531.4436     -87.294242      18.787515      18.776851      5.8701796    
        36  -1531.4441     -82.513789      18.787068      18.776535      5.8704075    
        37  -1531.4443     -81.685032      18.786867      18.776393      5.8705105    
        38  -1531.4443     -81.645941      18.786856      18.776385      5.8705163    
Loop time of 5.61363 on 1 procs for 38 steps with 192 atoms

1485.6% CPU use with 1 MPI tasks x 1 OpenMP threads

Minimization stats:
  Stopping criterion = energy tolerance
  Energy initial, next-to-last, final = 
     -1515.96938637628  -1531.44433389589  -1531.44434754222
  Force two-norm initial, final = 939.4985 6.5838898
  Force max component initial, final = 594.09391 3.1017355
  Final line search alpha, max atom move = 3.1475601e-07 9.7628989e-07
  Iterations, force evaluations = 38 72

MPI task timing breakdown:
Section |  min time  |  avg time  |  max time  |%varavg| %total
---------------------------------------------------------------
Pair    | 5.6035     | 5.6035     | 5.6035     |   0.0 | 99.82
Neigh   | 0.0014525  | 0.0014525  | 0.0014525  |   0.0 |  0.03
Comm    | 0.0017821  | 0.0017821  | 0.0017821  |   0.0 |  0.03
Output  | 0.0012273  | 0.0012273  | 0.0012273  |   0.0 |  0.02
Modify  | 0          | 0          | 0          |   0.0 |  0.00
Other   |            | 0.005695   |            |       |  0.10

Nlocal:            192 ave         192 max         192 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost:           2952 ave        2952 max        2952 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs:              0 ave           0 max           0 min
Histogram: 1 0 0 0 0 0 0 0 0 0
FullNghs:        53632 ave       53632 max       53632 min
Histogram: 1 0 0 0 0 0 0 0 0 0

Total # of neighbors = 53632
Ave neighs/atom = 279.33333
Neighbor list builds = 1
Dangerous builds = 0
#
#######################################################################
# (B) Final atomic positions refinement
 minimize 1e-8 1e-8 10000 100000
Generated 0 of 1 mixed pair_coeff terms from geometric mixing rule
WARNING: Energy due to 3 extra global DOFs will be included in minimizer energies
 (src/min.cpp:219)
Per MPI rank memory allocation (min/avg/max) = 4.301 | 4.301 | 4.301 Mbytes
   Step         PotEng         Press            Lx             Ly             Lz      
        38  -1531.4443     -81.645941      18.786856      18.776385      5.8705163    
        39  -1531.4444     -81.124708      18.78685       18.776381      5.8705192    
Loop time of 1.4188 on 1 procs for 1 steps with 192 atoms

1539.1% CPU use with 1 MPI tasks x 1 OpenMP threads

Minimization stats:
  Stopping criterion = energy tolerance
  Energy initial, next-to-last, final = 
     -1531.44434754222  -1531.44434754222  -1531.44435443502
  Force two-norm initial, final = 6.5903261 6.5899132
  Force max component initial, final = 3.0766263 3.0768585
  Final line search alpha, max atom move = 1.5870671e-07 4.883181e-07
  Iterations, force evaluations = 1 18

MPI task timing breakdown:
Section |  min time  |  avg time  |  max time  |%varavg| %total
---------------------------------------------------------------
Pair    | 1.4156     | 1.4156     | 1.4156     |   0.0 | 99.77
Neigh   | 0.0013218  | 0.0013218  | 0.0013218  |   0.0 |  0.09
Comm    | 0.00046373 | 0.00046373 | 0.00046373 |   0.0 |  0.03
Output  | 0          | 0          | 0          |   0.0 |  0.00
Modify  | 0          | 0          | 0          |   0.0 |  0.00
Other   |            | 0.001433   |            |       |  0.10

Nlocal:            192 ave         192 max         192 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost:           2952 ave        2952 max        2952 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs:              0 ave           0 max           0 min
Histogram: 1 0 0 0 0 0 0 0 0 0
FullNghs:        53632 ave       53632 max       53632 min
Histogram: 1 0 0 0 0 0 0 0 0 0

Total # of neighbors = 53632
Ave neighs/atom = 279.33333
Neighbor list builds = 1
Dangerous builds = 0

Total wall time: 0:00:08
2

Minimization is a numerical process affected by many factors, and generally speaking it’s impossible to get exact zero force or pressure. Your output shows that the final pressure is about 81bar, which is essentially same as zero in molecular simulations for a system of this size (the typical pressure fluctuations in MD would be much larger). Also you can note that the change of lattice parameters in the last steps are on the order of 0.0001A, which should be way smaller than the accuracy you need.

btw, not directly related to the problem, but

I also tried the approach taken by one of the official tutorials here

This really could only be determined by lammps developers, but I don’t think they are “official”.

1 Like
3

I had intuitively thought of this at first, but then I couldn’t account for the relatively quick pressure convergence for the much larger 5 x 5 x 5 supercell, which is nearly 4x the size. That system, according to your reasoning, should have much larger pressure fluctuations, correct?

As for the lattice parameters, I was concerned that it not changing despite being at nonzero pressure meant that it wasn’t really equilibriated per se, i.e not at the global energy minimum.

Intuitively, if I think about 80 bar pressure, it is quite high. That is a typical operating pressure for industrial hydraulic machinery. Is it reasonable to be approximating it to zero?

4

That system, according to your reasoning, should have much larger pressure fluctuations, correct?

No, indeed larger systems = smaller fluctuations. For detailed proof you can refer to your favorite textbook in statistical physics, but here is an intuitive argument: imagine splitting the 5x5x5 supercell into 125 unit cells. Each unit cell has a pressure, some higher than the average and some lower. The pressure of the whole system is some kind of average among all unit cells. Mathematics tell us that the average of many independent random variables has lower variance than each variable itself. So the pressure of the whole system has lower variance.

Intuitively, if I think about 80 bar pressure, it is quite high. That is a typical operating pressure for industrial hydraulic machinery. Is it reasonable to be approximating it to zero?

This is a common question asked by many people just entered the molecular simulation scene. In MD simulations, what matters is the geometry (i.e. distance between atoms), as it is what determines the forces, energies, etc. So generally you should look at the variation of lattice parameters, which already converged well in your case. It’s common to have quite large pressures “by common standards” in MD simulations for solids and liquids, and usually they are fine. On the other hand, if you check the phase diagram etc., this level of pressure rarely have significant effects on the properties of dense solid/liquids.

1 Like
5

It is next to impossible to make any definite statement without investing a lot of effort and re-check all of your work. So we can only speculate.

One possible explanation for difficulties to converge tightly, is that the potential hypersurface may be rather rugged. This commonly happens with tabulated potentials, which includes machine learning potentials. Withe the latter, it is also possible that you have not done sufficient training of your potential leading to inconsistent behavior.

No. Local pressure fluctuations are the largest. The larger the system, the more the fluctuations overlap and thus partially cancel and thus get smaller. This has been discussed here many times. Just search the forum.

You have a solid that is rather incompressible, so even the tiniest geometry changes can lead to large pressure fluctuations. You cannot apply macroscopic reasoning to atomic simulations.

1 Like
6

@initialize @akohlmey thank you. From all of this, I gather the most likely explanation for this is the incompressible nature of the material I am working with. In any case, I shall go back and re-check my MLIP parameters, and if possible, train it for longer.

I also observed that for my converged “zero pressure” system, while the overall pressure was zero, one of my two norm force components was very large, which resulted in the altered local symmetry. Here, I actually managed to obtain consistent bond lengths across octahedra, which is a good sign.