Dear Lammps users,

I have a question regarding the short and long range coulombic interactions. My question is why the short interactions (interactions inside the cuttoff; V=qi x qj x erfc(alpha rij)/rij ) is always calculated as positive value and long range interactions as negative values? I have wrote a simple code to calculate the coulombic interactions according to above formula and I am getting sometimes positive and sometimes negative numbers as the ensemble of pair charges (neg/neg , pos/pos, neg/pos) move in the system.

Could anyone please give a hint about this?

Many thanks

Dear Lammps users,

I have a question regarding the short and long range coulombic interactions. My question is why the short interactions (interactions inside the cuttoff; V=qi x qj x erfc(alpha rij)/rij ) is always calculated as positive value and long range interactions as negative values?

How do you conclude that?

I have wrote a simple code to calculate the coulombic interactions according to above formula and I am getting sometimes positive and sometimes negative numbers as the ensemble of pair charges (neg/neg , pos/pos, neg/pos) move in the system.

Could anyone please give a hint about this?

Most likely you made a mistake somewhere.

Axel

Deal Axel,

Thanks for the respond. I have written a simple routine for only the **short range** interactions (I copied the routine below). I am running this piece of code for a liquid (each molecule consists of 3 atoms of charges 0.25,-0.5, 0.25).

In my calculations the potential become positive for pair charges of 0.25/0.25 and -0.5/-0.5; and become negative for pair charges of 0.25/-0.5. These negative and positive numbers cancel out each other and I end up with very small number which fluctuate around zero. I ran the same system with lammps and the short range is a large positive number ( ~+1800 kcal/mol for all snapshots). I am really wondering what am I doing wrong here!!

I thank you in advance for your help.

uElect_shortRrange=0.0; //stores the electrostatic potential

for (int i=0; i< nmolec; i++) //zero the force

Force [i] = {0}; //zero the force

for (int ii=0; ii< nmolec -1; ii++) //nmolec is the number of molecules

{

for (int jj=ii+1; jj< nmolec ; jj++)

{

for (int iiSite=0 ; iiSite< natom; iiSite++) //natom is the number of atoms in the molecules (all molecules have the same number of atoms)

{

indexI=iiSite+ii* natom ; //find the index of atoms on molecule ii

for (int jjSite=0; jjSite< natom ; jjSite++)

{

indexJ=jjSite+jj* natom ; //find the index of atoms on molecule jj

dr.x=rx[indexI] - rx[indexJ];

dr.y=ry[indexI] - ry[indexJ];

dr.z=rz[indexI] - rz[indexJ];

dr.x-=fast_nearest(dr.x); //find the shortest distance in the periodic box

dr.y-=fast_nearest(dr.y); //find the shortest distance in the periodic box

dr.z-=fast_nearest(dr.z); //find the shortest distance in the periodic box

rijSqr = dr.x*dr.x + dr.y*dr.y + dr.z*drz;

rij = sqrt (rr);

if ( rij < rCut )

{

Ewald_force_mag = charge[indexI] * charge[indexJ] * { erfc (alpha * rij) / rijSqr /rij + 2. * alpha * exp (- alpha* alpha * rijSqr ) / sqrt(pi) / rijSqr };

Ewald_force.x= force_mag * dr.x; //force vector due to short range ewald

Ewald_force.x= force_mag * dr.y;

Ewald_force.x= force_mag * dr.z;

Force[ indexI ].x += Ewald_force.x;

Force[ indexI ].y += Ewald_force.x;

Force[ indexI ].z += Ewald_force.x;

Force[ indexJ ].x -= Ewald_force.x;

Force[ indexJ ].y -= Ewald_force.x;

Force[ indexJ ].z -= Ewald_force.x;

uElect_shortRrange += charge[indexI] * charge[indexJ] * erfc (alpha* rij)/ rij;

}

}

}

}

}

________________________________

Sorry, there was a little typo in the parameters

Arshia

From: Arshia Movahed <[email protected]...>

To: "[email protected]"

<[email protected]>

Sent: Monday, February 20, 2012 1:42 PM

Subject: Re: [lammps-users] short range interactions in Ewald sum

Deal Axel,

Thanks for the respond. I have written a simple routine for only the **short

range** interactions (I copied the routine below). I am running this piece

of code for a liquid (each molecule consists of 3 atoms of charges

0.25,-0.5, 0.25).

In my calculations the potential become positive for pair charges of

0.25/0.25 and -0.5/-0.5; and become negative for pair charges of 0.25/-0.5.

These negative and positive numbers cancel out each other and I end up with

very small number which fluctuate around zero. I ran the same system with

lammps and the short range is a large positive number ( ~+1800 kcal/mol for

all snapshots). I am really wondering what am I doing wrong here!!

where is the lammps input?

difficult to tell why there is a difference,

if i cannot see the difference.

axel.

Dear Axel,

I have attached my data file + input file. I am only comparing the *short range coulombic* part (the vdw part are not important for me).

Many Thanks

Arshia

PS: I apologize I had sent this to your email instead of lammps users

in.class2 (1.59 KB)

test_dme.lammps05 (31.2 KB)

ok. thanks.

now that i see your lammps input, the explanation is very simple:

your test code doesn't seem to handle exclusions,

i.e. skip over bonded interactions when computing

the non-bonded contributions.

if you change the lammps input by replacing:

special_bonds lj/coul 0.0 0.0 1.0

with

special_bonds lj/coul 1.0 1.0 1.0

you'll get a negative real space coulomb energy.

of course, this will be a bad model for a real

simulation, but that is beside the point here.

axel.

Dear Axel,

Thank you very much for your respond.

But the way that for-loops have been set, they would *not* count any interaction on the same molecule. And since my molecules have only 3 sites, the bond and angle terms (bonded) *are* in fact excluded (and I do not have any torsion of courseâ€¦)

And yes you are right about changing the special bond in lammps; I did check that.

I have become very confused that why I can not get the same results are lammps!!! It seems that I am making a *big* mistake since the numbers I get are very very different (1800 vs ~ 0 kcal/mol!).

Many thanks

Aeshia

And another thing is that, I get neg *and* positive values (fluctuating around zero since numbers cancel out each other)!

Thanks again

And another thing is that, I get neg *and* positive values (fluctuating

around zero since numbers cancel out each other)!

i don't have the time to debug your code,

but my suggestion is to first disregard the

short range screening and compare against

lj/class2/coul/cut

with a suitable cutoff and make sure that

that part is working properly and then

check what happens when you re-introduce

the screening.

axel.

Dear Axel,

Are the Windows source code (to run on series) available (I know that the executable file is available on the website)? Something that one can easily debug?

(debugging the linux version is difficult for me since I do have the root privilege for our cluster and hence building the code in debug mode would make lots of difficultiesâ€¦)

Thanks

Dear Axel,

Are the Windows source code (to run on series) available (I know that the

executable file is available on the website)? Something that one can easily

debug?

the sources for windows are the same as for unix.

(debugging the linux version is difficult for me since I do have the

root privilege for our cluster and hence building the code in debug mode

would make lots of difficulties...)

this is *complete* nonsense.

you don't need root to compile in debug mode

(or compile *anything*).

my experience with compiling on different platforms

is rather that compiling on windows is a major PITA.

axel.