Fix External and FIRE Minimization

Hi All,

I am new to LAMMPS and have run into an interaction between Fix External and FIRE minimization that I do not understand. If you have any insight into the behavior, I would very much appreciate the help. :slight_smile:

Following the examples in the Fix External documentation, I am using a Python callback function with this general form. (Here function(xyz) takes a set of xyz coordinates and returns global potential energy, atomic forces, and the global stress tensor. The Fix External instance is named “external potential”.)

 def callback(caller, ntimestep, nlocal, tag, xyz, fext, args = []):
    # evaluate function at coordinates xyz
    pe, forces, virial = function(xyz, *args)
  
    # update energy
    caller.fix_external_set_energy_global("external_potential", pe)
  
    # update forces
    fext = forces
  
    # update stresses
    caller.fix_external_set_virial_global("external_potential", virial)

I am attempting to use this callback to minimize structure using the following LAMMPS script. My test structure is DFT optimized, so add random displacement to the atoms prior to minimization.

units               metal
atom_style          atomic
dimension           3
boundary            p p p
box                 tilt large
read_data           {lammps_data_path}

comm_style          brick
comm_modify         cutoff 10.0

fix                 external_potential all external pf/callback 1 1
fix_modify          external_potential energy yes
fix_modify          external_potential virial yes

displace_atoms      all random 1e-3 1e-3 1e-3 13131313
min_style           fire
min_modify          integrator verlet tmax 6.0
minimize            0.0 1e-10 1000 10000

When run, the optimizer proceeds for a single step, during which the energy is unchanged and declares the force tolerance stopping criteria met. Curiously, when the callback function is asked to print the coordinates and forces, I find that LAMMPS is providing the same atomic coordinates for each call (and therefore receiving the same forces back in response). My assumption is that this is not the expected behavior?

Setting up fire style minimization ...
  Unit style    : metal
  Current step  : 0
WARNING: No pairwise cutoff or binsize set. Atom sorting therefore disabled. (src/atom.cpp:2142)
xyz:  [ 7.81302058  6.45639942 11.9702351 ]
forces:  [-1.3538678e-04  2.6862907e+00  1.0541548e-03]
Per MPI rank memory allocation (min/avg/max) = 2.724 | 2.724 | 2.724 Mbytes
   Step          Temp          E_pair         E_mol          TotEng         Press
         0   0              0              0             -778.56042     -2526719
xyz:  [ 7.81302058  6.45639942 11.9702351 ]
forces:  [-1.3538678e-04  2.6862907e+00  1.0541548e-03]
xyz:  [ 7.81302058  6.45639942 11.9702351 ]
forces:  [-1.3538678e-04  2.6862907e+00  1.0541548e-03]
xyz:  [ 7.81302058  6.45639942 11.9702351 ]
forces:  [-1.3538678e-04  2.6862907e+00  1.0541548e-03]
         1   0              0              0             -778.56042     -2526719
Loop time of 0.235064 on 1 procs for 1 steps with 110 atoms

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

Minimization stats:
  Stopping criterion = force tolerance
  Energy initial, next-to-last, final =
     -778.560424804688  -778.560424804688  -778.560424804688
  Force two-norm initial, final = 0 0
  Force max component initial, final = 0 0
  Final line search alpha, max atom move = 0 0
  Iterations, force evaluations = 1 2

Thanks again for taking the time to read through. I’ve spent a couple of weeks attempting to troubleshoot but have found myself stumped here. Any suggestions would be appreciated!

I am not sure that this will work as expected.

What happens, if you replace this with?

    for i in range(nlocal):
        fext[i][0] = forces[i][0]
        fext[i][1] = forces[i][1]
        fext[i][2] = forces[i][2]

Really appreciate the help! You were spot on - that has fixed the issue.