Recently I tried to develop a modified style of Coulombic force by editing both the pair_coul_cut.cpp and pair_coul.cut.h, and after this I successfully compiled and built the lmp executable with no error thrown out. For this new pair style I name it coul/test.
Then I tried to output the values trying to validate that my new style is working, I used the scripts below. However, I found that the set of values generated in this table.txt somehow is still as the shape of the original Coulombic potential.
So I doubt the pair_write command has used the new formula encoded in the .cpp file. I don’t know whether it’s I wrote the .cpp file wrong or I used the pair_write command in the wrong way.
The pair_write command uses calls to the Pair::single() function to
generate the table (with the charges being set to 1 and 1 if not given
explicitly in the command).
When thus to validate your implementation you need to update
Pair::single() as well and then I would suggest that you compare the
result of pair_write
to a system with 2 atoms where one of them is moved at constant
velocity with fix move and then you record energy and force via
thermodynamic output of fix print.
Thank you so much for your suggestion. When I try to implement my pair style in the Pair::single(), since there are seldom comments in the original source code, I have to guess what’s the physical correspondence to each variable, and it’s likely I guess them wrong.
In order to do it right, could you kindly inform me what are the physical quantities for these variables, either in the “single” or in the “compute” section of the source code (for coul/cut):
“Single” Section:
forcecoul (I guess it’s the potential?)
fforce
phicoul
“Compute” Section:
forcecoul
fpair
ecoul and eflag
And these seems to be the most important ones in the code from my naive observation. Also, generally speaking, I have been thinking about this question all along for the internal logic in LAMMPS about computing whatever forces: Do the source codes specify the formula of potential and force both and evaluate them respectively, or the codes just specify the formula for potentials and do some internal math of derivative to actually get the forces.
Really appreciate it if you can make some further explanation! Thank you so much for your patient instructions.
i don't really understand what there is to guess and why you need to
ask the questions that you are asking.
you have plenty of working examples where you *know* the potential and
how to compute the forces from them and that are *very* well
understood, e.g. the lj/cut or the coul/cut pair style.
just reading through the source code and matching it with the theory
should give you the answer you are looking for.
there is most certainly no automatic generation of the derivative to
get the forces. that would only work well for symbolic math and that
in turn is extremely slow compared to compiled software and speed and
efficiency is important in classical MD codes.
phi is often used for potential functions, so any variable with "phi"
in it is likely related to that. similarly for force.
class variables that are defined in the header or the header of the
parent class should be looked up in those headers, since those do have
comments explaining what they do, if it is not very obvious from
looking at the code itself.
please keep in mind two more things when reading the source code: 1)
that you are computing forces and energies for a spherical potential
based on the distance, but then also need project these distance based
forces on a 3-vector of x-, y-, and z-coordinates to get the forces on
the atoms in the lab frame. and 2) that for computational efficiency
the use of sqrt(), divisions and the power function are avoided (or if
it cannot be entirely avoided, reduced to the minimum) and invariant
prefactors, especially if they are global or per-type, precomputed.
the pow() function from the C standard math library is particularly
expensive since it assumes a floating point exponent and thus consists
of a costly exponential and a logarithm computation, while it can be
computed *much* faster for integer arguments (some c++ compilers have
specialized overloaded prototypes in the std:: namespace for that, but
you can also look at the MathSpecial::powint() function)
it may be simpler to study this with a simpler, standalone code.
please find below links to a lecture and source code that I used not
so long ago at an HPC workshop for a multi-day parallelization and
optimization project. it uses a simple, single atom type cutoff
lennard-jones MD code where you can see the discussion of the minimal,
naiive implementation first and how it can be optimized (it shows the
low-hanging fruit, not the maximum possible optimization) and
parallelized. in LAMMPS you don't need to worry about the latter, as
the (MPI) parallelization is abstracted from the compute kernels and
part of the neighbor list construction and time stepping framework.
For 2-dimension problem, in this case, we hope that if R> A, then A=R, where R is an atom-type variable and A is a common variable. The corresponding codes are as follows,
variable A equal 10
variable R atom sqrt((x-0)(x-0)+(y-0)(y-0)) # define a radius distance from point (0, 0)
if “v_R > $A” then “variable A equal V_R”
But I found that the above code cannot work. Any suggestions are very appreciated.
apart from being syntactically incorrect, the code below doesn't make
sense: the if statement is comparing a per atom property, but sets a
global property. so which atom's distance from the origin would count
as the global value? or what is it exactly you want to achieve?
Dear Axel,
Thanks for your reply.
In my code, I plan to heat the atoms to different temperatures (TR) according to the distance ® of each atom from a fixed point (0, 0). The position-dependent temperature TR is 1000*(cos(PIR/A)+1)/2 when the distance R is less than or equal to A. But when this distance R is greater than A , TR is 1000(cos(PI)+1)/2. So here we need a judgement statement of R and A.
for example,
variable A equal 10
variable R atom sqrt((x-0)*(x-0)+(y-0)*(y-0)) # define a radius distance from point (0, 0)
if " v_R < $A" then "variable TR atom 1000*(cos(PI*v_R/v_A)+1)/2 " else "variable TR atom 1000*(cos(PI)+1)/2 "
I know R is an atom variable, but A is a global variable. So the judgement is incorrect. But how can I achieve this purpose.
Dear Axel,
Thanks for your reply.
In my code, I plan to heat the atoms to different temperatures (TR) according to the distance (R) of each atom from a fixed point (0, 0). The position-dependent temperature TR is 1000*(cos(PI*R/A)+1)/2 when the distance R is less than or equal to A. But when this distance R is greater than A , TR is 1000*(cos(PI)+1)/2. So here we need a judgement statement of R and A.
this is a per-atom choice so you cannot implement this using a global
workflow command, it needs to be part of an atom style expression.
with those you can use the fact that logicals expand to either 1.0 or
0.0.
for example,
variable A equal 10
variable R atom sqrt((x-0)*(x-0)+(y-0)*(y-0)) # define a radius distance from point (0, 0)
so from here on you need to do something like this:
variable TA atom (v_R<v_A)*v_temp1+(v_R>=)*v_temp2
with v_temp1 and v_temp2 being the two different expressions to
compute the temperature. you compute them both, but due to multiplying
with the logicals they will be either multiplied with 1.0 or 0.0.
variable A atom 50 # variable A equal 50 ??? I had tried these two codes.
variable R atom sqrt((x-0)(x-0)+(y-0)(y-0)) #
variable TR1 atom 300+1000*(cos(3.1415926*v_R/v_A)+1)/2
variable TR2 atom 300
variable Tr atom (v_R<=v_A)*v_TR1+(v_R>v_A)*v_TR2
But I found that, both logicals of (v_R<=v_A) and (v_R>v_A) are always 0.0. So I doubt that these logical expressions of v_R and v_A are incorrect.
i cannot reproduce that. when i add the following lines to the melt
example in LAMMPS:
variable A atom 15.0
variable R atom sqrt((x-0)*(x-0)+(y-0)*(y-0)+(z-0)*(z-0))
variable T atom 1000*(cos(PI*v_R/v_A)+1)/2
variable TR atom 300+(v_R<=v_A)*v_T
dump 1 all custom 50 dump.lammpstrj id x y z v_R v_T v_TR
i get exactly the expected behavior, i.e. the variable TR computes to
1300 for the atom at 0.0 0.0 0.0 and then I get values between that
and 300.0 for atoms where the R variable evaluates to less or equal
than 15. and 300 for all above.
most likely you have a typo in your input deck somewhere.
