Re: [apbs-users] bug in coulomb?

[Andrew R. <dre...@gm...>]

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 <Nat...@pn...>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:dre...@gm...]
> *Sent:* Monday, February 24, 2014 11:22 PM
> *To:* apb...@li...
> *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