Dear all,
I have been working on simulating particle coagulation using pair_style hybrid/overlay with gran/hertz/history and born potential. In my simulation, I require different pair coefficients for different separation distance (d) ranges. For example, when d = 41.7 - 59.7: coefficients 1, d = 59.7 - 200: coefficients 2 and d = 200 - 700: coefficients 3. However, we can only assign an outer cutoff for each pair_coeff command in lammps. I tried to assign pair coeffs with an opposite sign to compensate for the excessive forces in each distance range. As a result, each range only has the potential as required. However, no agglomeration was observed using this method even though significant agglomeration should occur with the original potentials.
Does anybody have any idea of what I did wrong with hybrid/overlay and could you provide some suggestions?
Please find attached my input file. The lammps version I am using is 23Jun-2022.
Any help would be appreciated!
Thank you,
Shaelyn
in.in (2.0 KB)
I would plot an analytical form of your potential energy with gnuplot, something like:
# Borne potential definitions:
f1(x)=A1*exp((sigma1-x)/rho1) - C1/x**6 + D1/x**8
f2(x)=A2*exp((sigma2-x)/rho2) - C2/x**6 + D2/x**8
f3(x)=A3*exp((sigma3-x)/rho3) - C3/x**6 + D3/x**8
# Compound potential:
p(x)= x>0 && x < 41.7 ? 0 : x>=41.7 && x<59.7 ? f1(x) : x>=59.7 && x<200 ? f2(x) : x>=200 && x<700 ? f3(x) : 0
# Borne force definitions:
df1(x)= 0 # please differentiate df1/dx and write the analytical expression here
dp(x)= x>0 && x < 41.7 ? 0 : x>=41.7 && x<59.7 ? df1(x) : x>=59.7 && x<200 ? df2(x) : x>=200 && x<700 ? df3(x) : 0
set table "born.dat"
p [0:700] p(x), dp(x)
Please write the derivative to get the force. Then use the pair_style table to read your tabulated potential.
Hi Otello,
Thank you for your suggestions. However, table potential doesn’t work well for my system and that’s why I switched to born potential and gran pair style and used hybrid/overlay to combine them. Do you have any idea why pair hybrid/overlay doesn’t implement the calculation properly for this case?
Cheers,
Shaelyn
Why?
LAMMPS will do what you tell it to do. It knows nothing about the science that you want to implement, so if you don’t get the expected results the issue is likely that you don’t tell LAMMPS what you want to do. So first and foremost, you need to debug if you are getting the interaction that you want.
If you have a piece-wise defined potential using pair style table is the way to go. This is because usually you also have to employ some kind of switching function to avoid to have a discontinuous potential function as well as a discontinuous force, both can be a big problem an lead to unphysical results. In the tools/tabulate folder are some python classes and examples for building rather complex potential tables.
It may also be possible to define a piecewise potential using pair style lepton. You just define all segments as a sum, but then multiply each segment with like (r < “lower cutoff”)*(r > “upper cutoff”) to make its value 0.0 outside the expected range. Using Lepton is rather slow (but probly still competitive compared to looping over the same neighbors several times), but it allows to very quick experimentation and with the pair_write command you can also easily create a table from the final settings.