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------------------------------
##---------------INITIALIZATION-------------------------------
units metal
dimension 3
boundary p s s
atom_style molecular
##---------------ATOM DEFINITION------------------------------
read_data NP_A_DNT_16.data
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
##---------------minimize---------------------------------
min_style sd
minimize 1.0e-4 1.0e-6 100 1000
##---------------SETTINGS-------------------------------------
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
##---------------OUTPUT--------------------------------------
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
##---------------RELAXATION--------------------------------------
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.
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.