I would like to implement a simulation that utilizes tables for the pair potentials for both the short range coulombic interactions and the dispersion interactions. The tables I made for the coulumbic interactions are short range and I would like for it to be combined with the ewald long range calculations using kspace. I have not been able to get both sets of tables to work in the input file.
What I have tried:
hybrid overlay table linear 1500 {cutoff} coul/long{cutoff}
pair modify pair coul/long table linear 1500
However, I get errors saying illegal pair_style. How should I implement both tables at the same time? How can I let lammps know that one set of the tables is for coulombic together with kspace ewald, and the other set is for dispersion?
Also, for the coulombic tables, I made the potentials simply from (q*q/r)*damping_factor, is this adequate for it to fit together with kspace ewald?
I have a bit of confusion about the documentation regarding using tables, found in pair_table:
If your tabulated potential(s) are designed to be used as the short-range part of one of the long-range solvers specified by the kspace_style command, then you must use one or more of the optional keywords listed above for the pair_style command. These are ewald or pppm or msm or dispersion or tip4p. This is so LAMMPS can insure the short-range potential and long-range solver are compatible with each other, as it does for other short-range pair styles, such as pair_style lj/cut/coul/long. Note that it is up to you to insure the tabulated values for each pair of atom types has the correct functional form to be compatible with the matching long-range solver.
What does “compatible with the matching long-range solver” entail? How can I ensure that it matches with the long range solver? The table I made of the short range potential is simply the original coulombic interaction dampened at short distances.
Below is what I think might be how it works:
In creating the table, I assume I can use the functional form of the real space ewald solver: 0.5qq/rerfc(alpha*r), where alpha is the real/reciprocal space partition parameter (1/cutoff). From what I understand, I must match this partition parameter/cutoff by making the kspace_modify gewald value the same, thus matching the cutoffs. Is this correct and is there more matching I need to check to make sure it is compatible?
Additionally, if I would like to use tables for both coulombic and dispersion potentials, can I combine the energies in 1 table to be used in the realspace short range? I was unable to find information about making 2 sets of tables compatible in the documentation.
since you are obviously struggling to understand how tables work in LAMMPS, you should start by doing a tabulation for a functional form that exists analytically in LAMMPS. this way you can compare what you are doing with the tables to the corresponding results from the original analytical form. to start right away with your target is asking for trouble.
when tabulating coulomb, there are many issues to worry about, e.g. for tables to work you must have a different atom type for each atom that has a unique combination of atom type and charge. when using a coul/* pair style, this is not required, i.e. you can have the same atom type (i.e. the same lj/cut or other potential) with different atomic (partial) charges. with tables, those atoms have to have different types.
when starting with tables, it is often a big help to build the first set of tables using the pair_write command and then compare those with the output of the analytical form (see above) and the output of your program generating the tables
i would start first with a small test system without coulomb, then add coul/cut, then add coul/long and kspace
for tabulation of the coloumb, you can create a combined table with lj/cut (or equivalent) and coul/long (or equivalent) or create separate tables. in the second case, for hybrid/overlay, you have to specify pair style table twice.
there is little point in answering your questions unless you have worked your way through these suggestions and when you are done, it is likely, that you don’t need to ask those questions again.
Thank you! I have followed your suggestions to better understand how to use tables and have a few questions regarding that:
I first tested coul/cut with pair write to generate a table for that potential. I then ran the same simulation using pair_style table (with the tabulated potential of coul/cut), expecting to produce the same potential as coul/cut. The two simulation did produce the same energy, however the coul/cut approach allowed lammps to recognize that the energy was categorized as E_coul in the thermo output, while using pair_style table the energy was calculated as E_vdw. Is there a way to tell lammps that the tabulated potential using pair_table is coulombic energy rather than vdw?
I was able to match the analytical form of coul/cut with the pair write output of lammps. However, for the coul/long and kspace, what is the convergence parameter alpha (gewald in the source code?) that I can use in the analytical expression, following the functional form of the real space ewald pair wise calculation?
Thank you! I have followed your suggestions to better understand how to
use tables and have a few questions regarding that:
- I first tested coul/cut with pair write to generate a table for that
potential. I then ran the same simulation using pair_style table (with the
tabulated potential of coul/cut), expecting to produce the same potential
as coul/cut. The two simulation did produce the same energy, however the
coul/cut approach allowed lammps to recognize that the energy was
categorized as E_coul in the thermo output, while using pair_style table
the energy was calculated as E_vdw. Is there a way to tell lammps that the
tabulated potential using pair_table is coulombic energy rather than vdw?
no. the basic logic in LAMMPS is the following: if a pair-wise(!)
force/energy is computed based on the atom type it will be tallied into
E_vdw, if it is computed based on the individual atom's charge, it will be
tallied into E_coul. for tables - as explained previously - you *have* to
tie the charge to the atom type, since the charge value is included into
the table. thus the logical way to handle this is to tally the energy into
E_vdw.
- I was able to match the analytical form of coul/cut with the pair write
output of lammps. However, for the coul/long and kspace, what is the
convergence parameter alpha (gewald in the source code?) that I can use in
the analytical expression, following the functional form of the real space
ewald pair wise calculation?
it is automatically determined and should be printed to the screen when
you run a coul/long style. or you can set it explicitly yourself. for
details, please consult the manual.
I have implemented the coul/long tables and correctly matched analytically tabulated tables with the pair write tables produced from simulation as an exercise.
I am now trying to implement another functional form using tables. The functional form is the normal coulombic form with an erfc-style damping function and Buckingham-like dispersion contribution. I tabulated the table analytically using the functional form and used pair table to read the tables in the input file. I then used pair_write to check that the simulation is correctly reading the tabulated tables. However, there is discrepancy at short range (<0.5 angstrom) where the pair_write outputs the force converging to 0 (see attached graph) while my analytically tabulated tables go to a very high energy at small distances.
I found another questions in the mailing list archives that might be relevent:
where you mention that there is something regarding matching the cutoff that is used in screening with real space with the cutoff in screening. Can you clarify on what this means? The cutoff of my tabulated potential is 14 A, and inner cutoff of 0.1 A. The conversion parameter in the real space coulombic potential (alpha) I used is 0.1, which is from the gewald output for a coul/long simulation I ran with this model (this is perhaps not the exact alpha I should use, but I thought it would be a close estimate).
Could the discrepancies of my input tables and pair write tables be a result of cutoffs? What else can I check to figure out why the pair write is not outputting the tables I input?
I have implemented the coul/long tables and correctly matched analytically
tabulated tables with the pair write tables produced from simulation as an
exercise.
I am now trying to implement another functional form using tables. The
functional form is the normal coulombic form with an erfc-style damping
function and Buckingham-like dispersion contribution. I tabulated the
table analytically using the functional form and used pair table to read
the tables in the input file. I then used pair_write to check that the
simulation is correctly reading the tabulated tables. However, there is
discrepancy at short range (<0.5 angstrom) where the pair_write outputs the
force converging to 0 (see attached graph) while my analytically tabulated
tables go to a very high energy at small distances.
I found another questions in the mailing list archives that might be
relevent:
where you mention that there is something regarding matching the cutoff
that is used in screening with real space with the cutoff in screening.
Can you clarify on what this means? The cutoff of my tabulated potential
is 14 A, and inner cutoff of 0.1 A. The conversion parameter in the real
space coulombic potential (alpha) I used is 0.1, which is from the gewald
output for a coul/long simulation I ran with this model (this is perhaps
not the exact alpha I should use, but I thought it would be a close
estimate).
Could the discrepancies of my input tables and pair write tables be a
result of cutoffs?
unlikely. from your graph it looks like your analytical table has a
discontinuity at the lower cutoff, where the potential energy drops to
zero.
the output from pair_write looks like it was the result of fitting splines
to this (badly formed?) input. getting things right at the lower cutoff is
a known difficulty of using tables.
your data points for force and energy have to be smoooooth.
I hope to use a table for both the Coulomb and vdw potentials. You mentioned in an earlier email that I would need to use hybrid overlay to specify pair_style table twice. Can you clarify on how this is done? I used
pair_style hybrid overlay table linear 1400 ewald table linear 1400
but lammps doesn’t recognize this pair_style and I am also not sure how to tell lammps which table is which. How should the 2 separate tables (coulombic and vdw) be implemented?
Alternatively, can I combine the two tables by simply adding the potentials together to form 1 table? How would this work with the long range ewald summation?
I hope to use a table for both the Coulomb and vdw potentials. You
mentioned in an earlier email that I would need to use hybrid overlay to
specify pair_style table twice. Can you clarify on how this is done? I
used
pair_style hybrid overlay table linear 1400 ewald table linear 1400
it is "hybrid/overlay" not "hybrid overlay"
but lammps doesn't recognize this pair_style and I am also not sure how to
tell lammps which table is which. How should the 2 separate tables
(coulombic and vdw) be implemented?
