Fix bond/break and bond/create not working as expected

I need to have LAMMPS determine and display the number of bonds and its change every timestep. That isn’t possible if you use pairwise potentials, so I used bonds of “style zero”. I also used “fix bond/break and fix bond/create” to allow the change in the number of bonds according to the distance between the atoms and display it as I need. However, the bonds were only breaking, and I received a warning message saying “Fix bond/create is used multiple times or with fix bond/break - may not work as expected”. From this message, I supposed that the problem is from combining “fix bond/break and fix bond/create” in the same script, so I started looking for another way to get around this issue and reach my aims. Then I found “fix bond/react” which I don’t know if it is my best way for that. Is there any suggestions to fix this issue?

Please explain why you need to do this analysis during the run and at every timestep?

I want to determine the critical temperature at which the system becomes unstable. I consider the system unstable if the broken bonds reach a specific number after relaxation. Therefore, I added the “fix bond/break” command at first without the bond/create. However, I realized that this is not realistic because the bonds break and also form. So here is what I want from LAMMPS; start from an initial temperature, determine if the system is stable from the number of net broken bonds after each relaxation (iteration) and decide the next temperature accordingly until it reaches the critical temperature using the bisection numerical method.

it’s not very important to determine it every timestep, but it would be preferable. what really matters is to have it every relaxation.

I don’t see how you can do this from within LAMMPS automatically unless you program your own compute/fix to do the analysis report back the quantity you need to decide what to do next.

However, it should be rather straightforward to do this by doing a regular simulation where you also output the pairwise distances from compute pair/local and then analyze that with a custom script/program and prepare the next input based on the results from that.

The problem with fix bond/break and fix bond/create is that they cannot exchange information between what bonds were deleted or added by each of them during the same timestep. If you absolutely feel you must do the analysis during the run, I would consider compute coord/atom. It will count the number of neighboring atoms within a given cutoff. You would then compute those separately for C-C and C-H pairs, sum those per atom and over all atoms and then divide by 2.

I’ve read about the commands you suggested

I don’t think this will produce correct results since the atoms sometimes make bonds through the periodic boundery walls of the simulation box.

I believe that “compute pair/local” will be a very good and suitable option. It produces a vector which size is the number of the current bonds created by the pairwise interaction but leaving me confused on how to output the size value to go on with my analysis.

I don’t know if I came up with a good way to do so which is as follows; using the result of dividing “compute reduce sum” over “compute reduce ave”, but do you think this will produce a realistic and correct result?
also, do you know any other straightforward way to get the size of vectors in lammps?

Where does it say that compute coord/atom does not count neighbors across periodic boundaries?
It uses the same kind of neighbor list that is also used to compute pairwise interactions and those work fine across periodic boundaries.

As I already mentioned, the data from compute pair/local must be post-processed. This is because you cannot provide a cutoff distance. Instead it will use the neighborlist cutoff which is usually the largest pairwise cutoff plus the neighborlist skin distance. In post-processing the output form ‘dump local’ however, it is straightforward to determine the number of pair distances that are smaller than the desired cutoff value.

Perhaps the following is more to your liking:

variable d internal 0.0
variable btoolong equal v_d>2.0
compute bid all property/local batom1 batom2 btype
compute blen all bond/local dist v_btoolong set dist d
compute broken all reduce sum c_blen[2]
compute average all reduce ave c_blen[1]
thermo_style   custom step temp  etotal c_broken c_average
dump           2 all local  1000 NP_A_DNT_16_900.local c_bid[1] c_bid[2] c_bid[3] c_blen[1] c_blen[2]

This uses some special feature in the bond/local compute where you can post-process the data during data collection. While looping over the bonds, the distance is set to the variable d and then the value of btoolong is computed, which should be 1 for distances over 2 Angstrom and 0 for less. Then the local data from the compute can be subject to a compute reduce where you can sum the value of v_btoolong for all bonds and thus get the number of broken bonds. For demonstration purposes, I’ve also computed the average bond length (which only requires the “dist” property from compute bond/local). For verification, I’ve added a dump local command to output the identity of each bond and its distance as well as the current value for whether the bond is considered broken or not. The reduced values can also be output with the thermo_style command, and you’ll get something like this:

Step Temp TotEng c_broken c_average 
       0          900   -6594.8298            0    1.5233928 
     100    671.86567   -6583.4243            0    1.5495913 
     200    603.10702   -6564.8837            3    1.5645938 
     300    918.40214   -6530.7502           14    1.5572415 
     400     933.1734   -6526.7218           18    1.5491849 
     500    974.60762   -6537.7878           21    1.5612784 
     600    915.99415   -6551.8258           22    1.5706482 
     700    1033.2772   -6563.8236           29    1.5621487 
     800    1163.7881   -6588.4317           35    1.5636539 
     900    1209.9149   -6625.1785           39    1.5709774 
    1000    1165.2005   -6658.8034           47    1.5771916 
    1100    1127.0975   -6694.7238           49    1.5711762 
    1200    986.11058    -6722.003           51    1.5678127 
    1300    1005.8659   -6745.4795           54    1.5712717 
    1400    917.23451   -6763.5673           55    1.5748047 
    1500    924.46765   -6775.2171           56    1.5716404 
    1600    873.76411   -6783.2541           57    1.5705439 
    1700    899.33773   -6787.0693           57    1.5705623 
    1800    852.58721   -6787.8506           57    1.5729025 
    1900    926.57096   -6786.9968           59    1.5707542 
    2000    909.10739   -6792.1777           60    1.5686871 
    2100    893.17141   -6798.0091           61     1.570265 
    2200    897.84524   -6802.4336           61    1.5728448 
    2300    906.12195   -6803.9804           61    1.5733236 
    2400      929.719   -6805.1611           63    1.5699715 
    2500    981.75372   -6814.9155           68     1.576864 
    2600    1001.3266   -6826.1576           71    1.5759284 
    2700    1008.0125   -6837.2365           74    1.5753728 
    2800    1011.8262   -6846.2121           80    1.5782036 
    2900    947.92791   -6862.8091           82    1.5793407 
    3000    932.35875   -6876.6222           82    1.5783115 
    3100    926.32907   -6887.5912           84    1.5822947 
    3200    925.39344    -6898.616           87    1.5815561 
    3300    935.57056   -6908.5262           89    1.5816934 
    3400    961.43506   -6913.2645           90    1.5790956 
    3500    960.02943   -6919.7051           97    1.5885813 
    3600    999.78227   -6934.7299          101    1.5848615 
    3700    917.72128   -6947.2871          101    1.5829186 
    3800    910.64529   -6956.4815          101    1.5852419 
    3900    904.28228   -6963.2805          101    1.5871066 
    4000    902.05898   -6967.0908          102    1.5854285 
    4100    867.10468   -6968.1572          101    1.5845553 
    4200    896.55883   -6965.7105          101    1.5832532 
    4300    975.70692   -6958.9841          103    1.5842648 
    4400    946.56979   -6961.2152          103    1.5861205 
    4500    907.46255   -6967.9458          101    1.5835836 
    4600    933.96367   -6970.4453          105    1.5851978 
    4700    901.97801   -6975.7966          105     1.586907 
    4800    919.26366   -6980.5372          106    1.5855969 
    4900    908.22643   -6983.6344          105    1.5871658 
    5000    877.04143    -6984.716          105    1.5873795 
    5100    955.65277   -6983.3745          109    1.5885785 

hello Alex

I also used this command as follows

##---------------ATOM DEFINITION------------------------------


units          metal
dimension 	    3 
boundary       p s s
atom_style 	molecular

##---------------ATOM DEFINITION------------------------------

group G1 type 1
group G2 type 2
group G3 type 3
group G4 type 4

##---------------FORCE FIELDS---------------------------------
bond_style harmonic
bond_coeff * 0.0 1.0

pair_style 	airebo 3.0
pair_coeff     * * CH.airebo C C C C

special_bonds lj/coul 1.0 1.0 1.0

##---------------neighbor list---------------------------------

neighbor       2.0 bin
neigh_modify   every 1


min_style      sd
minimize       1.0e-4 1.0e-6 100 1000


timestep       0.001
variable       stemp equal 200
velocity       all create ${stemp} 12345 units box dist gaussian 
fix            1 all npt temp ${stemp} ${stemp} 0.1 x 0 0 0.1 drag 0.1 
fix            2 all momentum 1 linear 1 1 1 angular rescale

compute 11 G1 coord/atom cutoff 2.0 1
compute 12 G1 coord/atom cutoff 2.0 2
compute 13 G1 coord/atom cutoff 2.0 3
compute 14 G1 coord/atom cutoff 2.0 4

compute 22 G2 coord/atom cutoff 2.0 2
compute 23 G2 coord/atom cutoff 2.0 3
compute 24 G2 coord/atom cutoff 2.0 4

compute 33 G3 coord/atom cutoff 2.0 3
compute 34 G3 coord/atom cutoff 2.0 4

compute 44 G4 coord/atom cutoff 2.0 4

compute inBonds all reduce sum c_12
compute midBonds all reduce sum c_23
compute outBonds all reduce sum c_34

compute 11bonds all reduce sum c_11
compute 22bonds all reduce sum c_22
compute 33bonds all reduce sum c_33
compute 44bonds all reduce sum c_44

compute 14bonds all reduce sum c_14
compute 13bonds all reduce sum c_13
compute 24bonds all reduce sum c_24

variable inBrokenBonds equal 864-c_inBonds
variable midBrokenBonds equal 288-c_midBonds
variable outBrokenBonds equal 864-c_outBonds

variable totalBrokenBonds equal v_inBrokenBonds+v_midBrokenBonds+v_outBrokenBonds
#variable recreatedBonds equal c_11bonds*0.5+c_22bonds*0.5+c_33bonds*0.5+c_44bonds*0.5 +c_14bonds+c_13bonds+c_24bonds 

variable grossBrokenBonds equal ${totalBrokenBonds}-${recreatedBonds}
compute            pePerAtom all pe/atom 
compute            Peall all pe 

thermo_style   custom step temp etotal v_grossBrokenBonds v_inBrokenBonds v_midBrokenBonds v_outBrokenBonds &
v_recreatedBonds c_11bonds c_22bonds c_33bonds c_44bonds c_14bonds c_13bonds c_24bonds
dump           1 all custom 1000 NP_A_DNT_16_1C id type x y z c_pePerAtom

thermo         200
thermo_modify   lost warn


run            50000

, but I receicved this error:

ERROR: Variable inBrokenBonds: Compute used in variable between runs is not current (src/variable.cpp:1378)
Last command: variable grossBrokenBonds equal ${totalBrokenBonds}-${recreatedBonds}

which indicates that the compute (I don’t know which one exactly, but I suppose it’s the compute coord/atom) is not invoked during the run. I went through the mailing list archives to look for similar issues in which you answered them that the compute command must be consumed during the run to invoke its value before it is used in Variables formulas. A solution that you offered to someone is to use run 0 post no. I used this command but it didn’t work.

The problem you have here is a conceptual one and is cause by this line:

You (still) don’t seem to be able to grasp the difference between variable expansion with $ and a variable reference with v_ which is crucial for this.

There are other conceptual issues with that input (like defining multiple compute coord/atom and compute reduce even though there is no need and thus causing a massive negative performance impact).

Have you tried out my suggestion based on compute bond/local?

That should be much more efficient than looking at coordination numbers since it only considers the pairs of atoms for which you have defined bonds and does not need to loop over all neighbors of atoms within the cutoff determined by the force field.

I have, and it’s truely an efficient way. Also, I learned a powerful way to count sorted bond lengths, which is a thing that I have to thank you about. However, I had to change the model; the atoms at both ends of the nanotube were dangling causing wrong description of the system I need. I should study a section of the tube so both ends of my model had to have bonds through the periodic walls, and I didn’t know how to define bonds in the data file through these walls to go on with my simulation. So I felt compute coord/atom to be better since it doesn’t need bonds definition.
Additionally, I noticed some atoms at the edge formed bonds with its same type atoms as a result of the AIREBO pairing effect, but they weren’t counted on the bond/local since they have no bonds defined between them.

You can define bonds between any two atom IDs; with periodic boundaries, LAMMPS will pick the closest image when compiling the bond list (same as it does with neighbors). You will get a warning about inconsistent image flags, but that can be ignored in this specific case. However, since in your model the bonds are simply for accounting, you also have to have the correct box length, so that also the non-bonded interactions are reaching properly across the boundaries.

1 Like

Dear Alex

I want to fix this issue because whenever I try to develop my input file’s capabilities, I encounter many errors and it becomes on the list of potential bags cauzing those errors. I tried to reduce the unnecessary repetition of such commands by rewriting it this way:

compute            22 G2 coord/atom cutoff 2.0 2 3 4
compute            33 G3 coord/atom cutoff 2.0 3 4
compute            44 G4 coord/atom cutoff 2.0 4 

compute            inBonds all reduce sum c_11[2]
compute            midBonds all reduce sum c_22[2] 
compute            outBonds all reduce sum c_33[2]

compute            11Bonds all reduce  sum c_11[1]
compute            22Bonds all reduce sum c_22[1]
compute            33Bonds all reduce  sum c_33[1]
compute            44Bonds all reduce sum c_44[1]

compute            14Bonds all reduce sum c_11[4]
compute            13Bonds all reduce sum c_11[3]
compute            24Bonds all reduce sum c_22[3]

I know it could be reduced even more, but I tried to do it gradually in order to make it easier spotting bags and errors. When I ran this, I received this error message:

ERROR: Compute reduce compute does not calculate a per-atom array
Last command: compute            44Bonds all reduce sum c_44[1]

I checked the mannual of lammps LAMMPS-compute_reduce. It states this compute does calculate a per-atom array provided with examples!!!
Could you please help me to find what I missed?

Also, I have some questions related to this topic but discussicng about other commands. However, I don’t know if it would be better to post a new topic or just mention it here since it’s an extension of this discussion.

The issue is not with compute reduce (which BTW does NOT compute a per-atom array but rather a global scalar or vector depending on input. the whole point about reducing like you do is to sum over all atoms) but with the input to the reduce command, so the output of the compute coord/atom. That can compute one or more per atom coordination numbers and thus return either a per-atom vector or per-atom array, respectively.

Please note the syntax of that command very carefully because your input is highly redundant. The first compute is computing all per-atom properties also included in the second and third (c_22[2] is the same as c_33[1] and c_22[3] is the same as c_33[2]). Similarly the second also includes the result of the third (c_33[2] is the same as c_44), which is only a per-atom vector instead of an array and thus causes the error message. thus referring to c_44[1] is incorrect (as indicated) but only c_44 would be valid.

Hello Alex
I need to calculate the average bond lengths for unbroken bonds e.i. average with excluding the bonds with length greater than the cutoff 2A. Therefore, I used the way you suggested for me in this toppic

variable d internal 0.0
variable btoolong equal v_d>2.0
compute Bid all property/local batom1 batom2 btype
compute Blen all bond/local dist v_btoolong set dist d
compute unBrokenBondsLength all reduce min c_Bid[1] c_Bid[2] c_Bid[3] c_Blen[1] c_Blen[2] &
replace 1 5 replace 2 5 replace 3 5 replace 4 5

compute BondLengthAve all reduce ave c_unBrokenBondsLength[4]

However, I received this error

ERROR: Compute reduce compute calculates global values (src/compute_reduce.cpp:214)
Last command: compute BondLengthAve all reduce ave c_unBrokenBondsLength[4]

does this mean I can’t use this method to apply the statistical analysis selectively for unborken bonds?
Do you think this can be addressed with some adjustments on the code? or do you suggest another way to achieve my goal.

What you are doing makes syntactically and logically no sense.
You need to get better at reading documentation and logical thinking. At this point I have already spent too much time on this with too little constructive input from you and you have exhausted my willingness to do your thinking and your work for you.