About adding new potential

LAMMPS LAMMPS Mailing List Mirror
1

Hi users,

I am trying to add one new potential for my physics problem in plasmas.
Just to start with (I have read the modification section in manual), I choose one coul/cut pair potential and just changed the name of potential
as coul/cut/on and similarly changed the CLASS in .cpp and .h files.
When I recompile lammps, it goes fine and then even I can run a test input file with this new potential (Definitely results are same as in coul/cut as no change made).

Further, I just made a small change (say added a factor while calculating fpair) and again after recompiling, the new potential is giving some result, different from coul/cut.

But when I pair_write this new potential and compare it with analytic
coulomb potential, both are same.

I checked it few many times …but unable to get error !!

What could be the probable error !!

Thanks

2

Hi users,

I am trying to add one new potential for my physics problem in plasmas.
Just to start with (I have read the modification section in manual), I
choose one coul/cut pair potential and just changed the name of potential
as coul/cut/on and similarly changed the CLASS in .cpp and .h files.
When I recompile lammps, it goes fine and then even I can run a test input
file with this new potential (Definitely results are same as in coul/cut as
no change made).

Further, I just made a small change (say added a factor while calculating
fpair) and again after recompiling, the new potential is giving some result,
different from coul/cut.

But when I pair_write this new potential and compare it with analytic
coulomb potential, both are same.

I checked it few many times ..but unable to get error !!

What could be the probable error !!

you changed the PairXXX::compute() method (which processes neighbor
lists) but not the PairXXX::single() method (which operates on giving
atom types and (squared) distance). the pair_write command utilizes
the latter, hence the discrepancy you see.

axel.

3

Hi Axel,

Thanks for always being prompt ! I already checked PairXXX::single() which is using same rsq (dx^2 + dy^2 + dz^2)
which has been calculated in PairXXX::compute() and I have made changes only in this “rsq” variable to make required changes for my potential.

Anyway, I am checking again !!

Thanks

4

Further, I just made a small change (say added a factor while calculating
fpair) and again after recompiling, the new potential is giving some result,
different from coul/cut.

But when I pair_write this new potential and compare it with analytic
coulomb potential, both are same.

  A quick comment: Are you saying that BOTH energy AND the force have
not changed since you altered the code? (As opposed to checking only
the energy. The forces are stored in the 4th column of the file
created by pair_write.)

Keep in mind that the code which calculates the forces, and the code
which calculates the energies are typically in different parts of the
PairXXXX::compute() function, and can be totally incompatible with
each other. (For example the energy for pair_style coul/cut is
calculated in "pair_coul_cut.cpp" near line 112. [2015-1-21 version]
"if (eflag)"...)

This is not the only explanation for your problem, but it seems like a
likely one.

Andrew

5

Further, I just made a small change (say added a factor while calculating
fpair) and again after recompiling, the new potential is giving some result,
different from coul/cut.

But when I pair_write this new potential and compare it with analytic
coulomb potential, both are same.

  A quick comment: Are you saying that BOTH energy AND the force have
not changed since you altered the code? (As opposed to checking only
the energy. The forces are stored in the 4th column of the file
created by pair_write.)

Keep in mind that the code which calculates the forces, and the code
which calculates the energies are typically in different parts of the
PairXXXX::compute() function, and can be totally incompatible with
each other. (For example the energy for pair_style coul/cut is
calculated in "pair_coul_cut.cpp" near line 112. [2015-1-21 version]
"if (eflag)"...)

This is not the only explanation for your problem, but it seems like a
likely one.

sorry, but the pair_write command doesn't use the compute() method at all.

please check out the Pair::write_file() method in pair.cpp around
lines 1465 and 1570.

axel.

6

you changed the PairXXX::compute() method (which processes neighbor
lists) but not the PairXXX::single() method (which operates on giving
atom types and (squared) distance).

Oops, yes, I forgot about that. Thanks Axel.

   So, Sanat, if you change the way you calculate pair interactions,
you have to modify 4 places in the code to keep the code consistent
with itself (to calculate the energies & forces for both
PairXXX::compute, and PairXXX::single()).
It's pretty easy to forget one of them.
Cheers

Andrew

7

Thanks Axel and Andrew !!

I got the point hopefully !! Let me implement it and make checks and then I will post the outcome.

Regards

8

Hi,

I just followed suggestions from Axel and Andrew and made changes to PairXXX::compute() method (which processes neighbor lists) and to the PairXXX::single() method both (in my case… added same constant to “rsq” variable in both functions).
Now pair_write is giving proper results.

Now I want to add this constant through “Pair_coeff” in input file. I looked over many pair_style.cpp files for how they give Pair coeff and then followed to make one for mine case. As a example… again coul/cut case …just need to add a more coeff for new pair potential…

Only changing things in PairCoulCutOn::coeff (given following) looks not enough as the error comes while compiling as:
“error: ‘lcut’ was not declared in this scope”

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

void PairCoulCutOn::coeff(int narg, char **arg)
{
if (narg < 3 || narg > 4)
error->all(FLERR,“Incorrect args for pair coefficients”);
if (!allocated) allocate();

int ilo,ihi,jlo,jhi;
force->bounds(arg[0],atom->ntypes,ilo,ihi);
force->bounds(arg[1],atom->ntypes,jlo,jhi);

double lcut = force->numeric(FLERR,arg[2]);

double cut_one = cut_global;
if (narg == 4) cut_one = force->numeric(FLERR,arg[3]);

int count = 0;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo,i); j <= jhi; j++) {
cut[i][j] = cut_one;
scale[i][j] = 1.0;
setflag[i][j] = 1;
count++;
}
}

if (count == 0) error->all(FLERR,“Incorrect args for pair coefficients”);
}

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
I am not familiar quite much with c++ but trying to catch logic.

Resolution will resolve my problem of this new potential. Thanks in advance.

9

The error message printed out by the compiler will tell you which line
it thinks the error is located on. I don't think it occurs in this
excerpt of text you posted because the "lcut" variable only appears
where it is declared.

Your problem appears to be that you have declared the double lcut
variable as a local variable inside PairCoulCutOn::coeff(). That
means other functions functions (like PairCoulCutOn::compute()) do not
have access to that variable.
You want to add a member variable to the PairCoulCutOn class. To do that:

1) Edit the "pair_coul_cut_on.cpp" file and replace:

double cut = force->numeric(FLERR,arg[2]);

  with:

cut = force->numeric(FLERR,arg[2]);

10

Only changing things in PairCoulCutOn::coeff (given following) looks not
enough as the error comes while compiling as:
"error: ‘lcut’ was not declared in this scope"

The error message printed out by the compiler will tell you which line
it thinks the error is located on. I don't think it occurs in this
excerpt of text you posted because the "lcut" variable only appears
where it is declared.

Your problem appears to be that you have declared the double lcut
variable as a local variable inside PairCoulCutOn::coeff(). That
means other functions functions (like PairCoulCutOn::compute()) do not
have access to that variable.
You want to add a member variable to the PairCoulCutOn class. To do that:

1) Edit the "pair_coul_cut_on.cpp" file and replace:

double cut = force->numeric(FLERR,arg[2]);

  with:

cut = force->numeric(FLERR,arg[2]);

...continuing:

2) Declare that "lcut" is a data member of class PairCoulCutOn, by
editing the "pair_coul_cut_on.h" file and add

double lcut;

somewhere in the list of "public:" members.

Editing the LAMMPS code requires some basic knowledge object-oriented
programming. Take a look at these web pages before continuing:

http://www.cplusplus.com/doc/tutorial/classes/
http://www.learncpp.com/cpp-tutorial/89-class-code-and-header-files/

Andrew
(I hate it when my message gets truncated. Sorry about that. Please
do visit those web pages before posting other questions about C++.)