## [apbs-users] bug in coulomb?

 [apbs-users] bug in coulomb? From: Andrew Ritchie - 2014-02-25 07:22:04 Attachments: Message as HTML ```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 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 , 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; ialist); iatom++) { atom = Valist_getAtom(thee->alist, iatom); apos = Vatom_getPosition(atom); charge = Vatom_getCharge(atom); for (ipos=0; ipos 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 ```

 [apbs-users] bug in coulomb? From: Andrew Ritchie - 2014-02-25 07:22:04 Attachments: Message as HTML ```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 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 , 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; ialist); iatom++) { atom = Valist_getAtom(thee->alist, iatom); apos = Vatom_getPosition(atom); charge = Vatom_getCharge(atom); for (ipos=0; ipos 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 ```
 Re: [apbs-users] bug in coulomb? From: Baker, Nathan - 2014-02-26 01:28:45 Attachments: Message as HTML ```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@...] Sent: Monday, February 24, 2014 11:22 PM To: apbs-users@... 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 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 , 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; ialist); iatom++) { atom = Valist_getAtom(thee->alist, iatom); apos = Vatom_getPosition(atom); charge = Vatom_getCharge(atom); for (ipos=0; ipos 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 ```
 Re: [apbs-users] bug in coulomb? From: Andrew Ritchie - 2014-02-26 01:39:22 Attachments: Message as HTML ```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 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@...] > *Sent:* Monday, February 24, 2014 11:22 PM > *To:* apbs-users@... > *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 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 , 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 > 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; iatomalist); iatom++) { > > atom = Valist_getAtom(thee->alist, iatom); > > apos = Vatom_getPosition(atom); > > charge = Vatom_getCharge(atom); > > for (ipos=0; 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 ```
 Re: [apbs-users] bug in coulomb? From: Fligg, Keith - 2014-02-26 01:51:06 Attachments: Message as HTML ```Will do! - Keith On Feb 25, 2014, at 18:28, "Baker, Nathan" > 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@...] Sent: Monday, February 24, 2014 11:22 PM To: apbs-users@... 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 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 , 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; ialist); iatom++) { atom = Valist_getAtom(thee->alist, iatom); apos = Vatom_getPosition(atom); charge = Vatom_getCharge(atom); for (ipos=0; ipos 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 ```