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. 
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!