Nathan,

After thinking about this some more it occurred to me that I was interpreting the reported forces as the force acting on each atom when the actual intention may be (probably is?) to report the force each atom is exerting (Newton's third law), in which case the sign convention is correct.     

Andrew


On Tue, Feb 25, 2014 at 7:28 PM, Baker, Nathan <Nathan.Baker@pnnl.gov> wrote:

Hi Andrew –

 

Thanks for sending this.

 

Keith, can you please look into this and patch the code, if needed?

 

Thanks,

 

--

Nathan Baker

Laboratory Fellow and Technical Group Manager

Applied Statistics and Computational Modeling

Pacific Northwest National Laboratory

+1-509-375-3997

http://www.linkedin.com/in/nathanandrewbaker/

 

From: Andrew Ritchie [mailto:drew.w.ritchie@gmail.com]
Sent: Monday, February 24, 2014 11:22 PM
To: apbs-users@lists.sourceforge.net
Subject: [apbs-users] bug in coulomb?

 

If I put a +2 charge at 0,0,0 and a +1 test charge at 1,0,0 and use apbs/tools/bin/coulomb -f <pqr> then the reported force on the +2 charge is in the positive x-direction and the reported force on the +1 charge is in the negative x-direction, which are the opposite of the expected signs.  I'm running the version which came with APBS 1.4.1.

 

Looking through coulomb.c, forces are being calculated simply as F=qE where there is an additional scalar -0.5.  I believe the scalar -0.5 should just be 0.5.  I looked though the Vgreen_coulombD_direct function in vgreen.c to make sure it wasn't calculating the negative field and thus needed to multiply by negative to return the correct force, and my first thought looking through the code is that's exactly what was happening.  However, closer inspection led me to realize that the outer loop is the atom contributing the field and the inner loop is the field location, thus the vectors dx,dy,dz are negative and the -= sign makes sense without the -0.5 scalar in coulomb.c.

 

I would just like to verify that this is correct before I try to rely on these results.

 

I've copy/pasted the test .pqr text, output from coulomb -f <pqr>, and the relevant code-snippets below. 

==================================================

PQR:

ATOM      1 TCHG TCHG    1       1.000   0.000   0.000  1.00  0.000

ATOM      2  MG   MG        2       0.000   0.000   0.000  2.00  0.800

 

Coulomb:

/path/to/apbs/tools/bin/coulomb -f ion.pqr 

 

Providing per-atom forces...

Setting up atom list from ion.pqr.

Read 2 atoms

Setting up Green's function object.

Dielectric constant = 1 (vacuum)

Distances in Angstroms

Charges in electrons

Allocating space for solution...

Calculating...

        Atom 1:  x-force = -1.389354796845E+03 kJ/mol/A

        Atom 1:  y-force = -0.000000000000E+00 kJ/mol/A

        Atom 1:  z-force = -0.000000000000E+00 kJ/mol/A

        Atom 2:  x-force = 1.389354796845E+03 kJ/mol/A

        Atom 2:  y-force = -0.000000000000E+00 kJ/mol/A

        Atom 2:  z-force = -0.000000000000E+00 kJ/mol/A

 

 

-------------------------------------------------------

Total energy = 2.778709593689e+03 kJ/mol in vacuum.

Total x-force = 0.000000000000e+00 kJ/mol/A in vacuum.

Total y-force = 0.000000000000e+00 kJ/mol/A in vacuum.

Total z-force = 0.000000000000e+00 kJ/mol/A in vacuum.

 

==================== coulomb.c ===================

    for (i=0; i<Valist_getNumberAtoms(alist); i++) {

        fx[i] *= (-0.5*qp[i]*zmagic);

        fy[i] *= (-0.5*qp[i]*zmagic);

        fz[i] *= (-0.5*qp[i]*zmagic);

        ...

 

==================== vgreen.c ====================

 

    for (iatom=0; iatom<Valist_getNumberAtoms(thee->alist); iatom++) {

        atom = Valist_getAtom(thee->alist, iatom);

        apos = Vatom_getPosition(atom);

        charge = Vatom_getCharge(atom);

        for (ipos=0; ipos<npos; ipos++) {

            dx = apos[0] - x[ipos];

            dy = apos[1] - y[ipos];

            dz = apos[2] - z[ipos];

            dist2 = VSQR(dx) + VSQR(dy) + VSQR(dz);

            dist = VSQRT(dist2);

            if (dist > VSMALL) {

                idist3 = 1.0/(dist*dist2);

                gradx[ipos] -= (charge*dx*idist3);

                grady[ipos] -= (charge*dy*idist3);

                gradz[ipos] -= (charge*dz*idist3);

                pot[ipos] += (charge/dist);

            }

        }

    }

==================================================

 

All the best,

Andrew Ritchie




--
-Andrew Ritchie