Fix addforce command in minimization process

LAMMPS LAMMPS Beginners
1

Dear LAMMPS developers,

I am attempting to relax a two-dimensional system consisting of approximately 500 atoms using the minimize command with the min_style hftn. In this process, I am adding externally calculated per-atom forces and potential energy iteratively using the fix addforce command (in addition to the forces arising from a pairwise potential ‘rebo’ being used). After each iteration, I allow the minimization step to proceed, during which the atoms adjust their positions in response to the forces acting on them, with the total energy being minimized over several iterations.

Recently, I discovered a conceptual error in my theory, which led me to double the magnitude of the forces and potential energy I had been using for the same system. My expectation was that by doubling the magnitude of these quantities, the atoms would move more significantly than in the previous case, and the minimization would be achieved in fewer iterations. However, I observed that even after adding the doubled forces, the atoms did not move as expected.

I would appreciate your insights into why this might be the case. I am using LAMMPS version 2 Aug 2023 with the Python interface. Below is the LAMMPS input I am using.

Note: For each iteration, I update the forces to be added based on the information in the dump file created in each iteration.

for run in range(200):

    units metal                                        
    dimension 3                                        
    boundary p p p                                     
    atom_style full                                    

    read_data structure.mol

    newton on                                 

    pair_style hybrid rebo zero 5.0 full
    pair_coeff * * rebo CH.rebo  C C
    pair_coeff 1 2 zero 

    neighbor 2 bin                                      
    neigh_modify every 1 delay 0 check yes               

    timestep  0.1

    variable zero atom 0.0
    dump             1 all custom 1 dump1.minimization id  type x y z fx fy fz #🔎 to be consistent with the BLBLNew.mol file format

    # Below: The v_nFx v_nFy v_nFz are the peratom forces read as variables 
    # The v_naPE is the peratom potential energy read as a variable
    
    fix 3 all addforce v_nFx v_nFy v_nFz every 1 energy v_naPE 
    fix_modify 3 energy yes

    min_style  hftn
    minimize 1e-5 1e-5 1 10000
2

This is bogus. Pair style rebo is not a pairwise potential, but disabling pairs between atom type 1 and atom type 2 you are also affecting the forces for triples and quadruples. If you do not want any interactions between atom types 1 and atom types 2 but rebo among themselves you would have do do something like this:

 pair_style hybrid rebo rebo
    pair_coeff * * rebo 1 CH.rebo  C NULL
    pair_coeff * * rebo 2 CH.rebo  NULL C
3

Dear Axel,

Thank you for your explanation, and I apologize for not providing the context for the commands below. I work on layered materials where we have both interlayer and interlayer interactions.

pair_style hybrid rebo zero 5.0 full
pair_coeff * * rebo CH.rebo  C C     # intra-layer
pair_coeff 1 2 zero                  # interlayer

In my input, I needed to incorporate externally calculated forces and potential energy while minimizing the energy of the system, particularly for (interlayer) interactions between atoms of type 1 and 2. Atoms of Type 1 are in one layer and atoms of type 2 are in the other layer.

To achieve this, I utilized the ‘zero’ pair style for the interlayer interactions. This allowed me to generate neighbor lists (‘find_pair_neighlist(‘zero’)’) using this pair style, enabling me to calculate per-atom forces and potential energy based on these neighbor lists, which need to be added to respective atoms using the ‘fix addforce’ command. In my work, I follow an approximation model that relies solely on the closest single interlayer neighbor.

For interactions between atoms of type 1 interacting with other atoms of type 1 and type 2 interacting with other atoms of type 2, I opted for the ‘rebo’ pair style. By utilizing the ‘fix addforce’ command, I could add externally calculated per-atom forces to the existing per-atom forces generated by the ‘rebo’ pair style.

My main confusion arises when doubling the amplitude of these external forces and potential energy (prefactor = 2). Initially, when the prefactor = 1, the atoms move as expected. However, when I increase the prefactor beyond 1.4, the atoms suddenly stop moving altogether. I would have expected them to move more significantly.

These tests using gradually increasing prefactors make me think there is something “failing” at the level of the minimization scheme. Based on my limited understanding, as I am only doing one minimization step within LAMMPS using the htfn algorithm after updating the forces externally, the number of minimization steps, set to 1, sets here the outer loop number of iterations to 1.

The inner loop, which takes care of finding the direction in which to move the atoms, is still equal to the default 1000. If the atoms are not moving at all, as is the case when my prefactor > 1.4, I would maybe expect some lammps error saying that the linesearch alpha is zero. However, I am just getting the usual message that LAMMPS finishes because the number of maximum (outer loop) iterations has been reached. Because of my slightly unconventional way to use LAMMPS, i.e. doing a single LAMMPS outer loop iteration interfacing with my forces calculated and updated externally, before doing the next single LAMMPS outer loop iteration, etc… maybe there is something that I am missing or unaware of on how LAMMPS generates such warnings. Do you have any pointer?

Or do you have any other suggestions on what could cause the atoms not to move at all if the forces are too strong? And possible ways to resolve them?

4

You have not paid attention to what I have explained and your pair style and pair coeff setup remains bogus. If you need to add pair style zero to get a neighbor lists, you can add that to the my suggested pair style and pair coeff setup. But only with the way that I was showing you can you have only intralayer interactions modeled by rebo.

5

I don’t really understand what you are trying to do with the workflow you have designed. It make no sense to me starting from the point that you limit the minimizer to a single step.