From: Nathan A. Baker <baker@ch...>  20040224 22:46:05

Dear Peter  Thanks for your email. The problems you are encountering are due to grid artifacts. Specifically, you are attempting to use the multigrid solver to numerically evaluate the (very) singular Green's function. This causes several problems: First, the exact Green's function contribution to the forces and energies (the "selfforces" and "selfenergies") is infinite and should be removed for a meaningful measurement of the system. Second, numerical solvers have a very hard time resolving the rapidlychanging Green's functions near the delta function sources. The resulting solution is very sensitive to the grid parameters such as the placement of the delta functions and the grid spacing. The result of all this is the artifical net force on the system in the absence of any external field. As an example, try setting "cgcent" and "fgcent" to specific points (rather than mol 1) and you'll observe significant variation in the forces. This is all a symptom of inaccuratelycalculated, very large, selfforces that do not cancel correctly. PB solvers, such as APBS, are generally used to determine solvation energies and forces for molecules  not the Coulombic types of interactions you were attempting calculate. You should find very stable forces if you calculate the difference between the "pdie 1.0 sdie 1.0" reference calculation you were performing and a "pdie 1.0 sdie 80.0" solvated calculation. Thank you for your interest in APBS!  Nathan Baker Peter Fedichev <PeterOFedichev@...> (02242004 11:49:300500): >Dear Nathan Baker, > >first of all thanks for creating, supporting and sharing the APBS software. > >I am evaluating the APBS solver for possible applications in >electronic structure calculations and up to the moment could only go >as far as to running simple examples. The simplest test, of course, >is the Coulomb law. Below I provide a simple complex.prq file >containing only 3 charges: > >ATOM 1 I ION 1 3.000 0.000 0.000 1.000 1.000 >ATOM 2 I ION 1 2.000 0.000 0.000 1.000 1.000 >ATOM 3 I ION 1 0.000 2.000 0.000 0.100 1.000 > >I am setting up the simplest vacuum calculation (epsilons=1). The >result of the calculations is checked first by calculating the total >force acting on the "molecule" expecting that the sum of all forces >in a closed system vanishes. > >I have tried both mgmanual and mgauto, my current settings are: > >elec > mgauto # Use the multigrid method > dime 65 65 65 # Grid dimensions (c 2^(l+1) +1) > # (in x, y, z) > nlev 4 # Number of levels in multilevel > # hierarchy (usually 'l' from above > # definition) > #grid 0.21 0.21 0.21 # Grid spacings (A) > glen 50. 50. 50. # You can EITHER specify the grid > cglen 300.00 300.00 300.00 << enormous values, see below > fglen 30. 30. 30. > # lengths OR the grid spacing > gcent mol 1 # This can either be 'mol n' to center > cgcent mol 1 > fgcent mol 1 > # on a specific molecule or (x, y, z) > # coordinates > mol 1 # Which molecule (1, 2, ...) > lpbe # lpbe/npbe = linear/nonlinear PBE > bcfl 1 # Boundary condition flag > # 0 => Zero > # 1 => Single DH sphere > # 2 => Multiple DH spheres > # 4 => Focusing > #ion 1 0.000 2.0 # Counterion declaration: > #ion 1 0.000 2.0 # ion <charge> <conc (M)> <radius> > pdie 1. # Solute dielectric > sdie 1. # Solvent dielectric > chgm 1 # Charge disc method (linear) > srfm 2 # Surface calculation method > # 0 => Mol surface for epsilon; > # inflated VdW for kappa; no > # smoothing > # 1 => As 0 with harmoic average > # smoothing > # 2 => Cubic spline > srad 1.4 # Solvent radius > swin 0.3 # Surface cubic spline window > temp 298.15 # Temperature > gamma 0.105 # Surface tension (in kJ/mol/A^2) > calcenergy 1 # Energy I/O (to stdout) > # 0 => don't write out energy > # 1 => write out total energy > # 2 => write out total energy and all > # components > calcforce 1 # Atomic forces I/O (to stdout) > > >*) enormous values for cglen parameter is my attempt to move the >boundaries as further away as i could. I have tried different values >of all *len parameters, the problem does not go away > >Now the results of the calculation: > > Global net fixed charge force = (2.514477272749E+00, 8.933277299877E+02, 4.760514198420E05) kJ/mol/A > >Which are quite big forces.... If I leave only the first two lines of my complex.pqr I get > > Global net fixed charge force = (2.551021300247E05, 4.465403497053E05, 4.612637773136E05) kJ/mol/A > >which is somewhat better (or maybe even exact, given the fact that >the whole system is not closed and boundaries could lead to the >overall nonvanishing contribution). Further calculation of the >forces acting on individual atoms does show that the forces are >calculated right. > >I am sorry for the long introduction, now the question: Why setting >up another (small) charge in a three charge system leads to such >dramatic change of the forces? Is it my fault? I have seen (maybe) >related issues in the APBS mailing list achieve in the past. I know >the suggested solution: calculate with epsilon=(say)80, then with >epsilon=1, subtract the two and add the result from the coulomb >binary. Unfortunately this did not help (at least in this example). >Having some (limited) experience with other solvers (numerical >recipes multigrid or FFTbased solutions for spatially homogeneous >systems to name a few) I do not understand why one should expect such >strong "selfforces" to be subtracted. > >That's been a lengthy complain, but I tried to provide as much info >as I could. Could you please suggest a direction to proceed with the >program? I am sure further clarification on this mailing list can be >of valuable help for other people intending to use the program in >their studies. > >Many thanks in advance >Peter > > End of message from Peter Fedichev .  Nathan A. Baker, Assistant Professor Washington University in St. Louis School of Medicine Dept. of Biochemistry and Molecular Biophysics Center for Computational Biology 700 S. Euclid Ave., Campus Box 8036, St. Louis, MO 63110 Phone: (314) 3622040, Fax: (314) 3620234 URL: http://www.biochem.wustl.edu/~baker PGP key: http://cholla.wustl.edu/~baker/pubkey.asc 
From: PeterOF<edichev@ne...>  20040225 18:09:43
Attachments:
apbs.in

Dear Nathan, thank you very much for your mail. I have looked through it carefully and still would like to ask you for further clarification. I do understand that any numerical solver is going to have selfenergies and selfforces problems. This is why I am aiming at the forces first rather than for the energies. I also know that pointlike charges are fundamentally problematic for multigrid methods. With all this in mind I did try to do what you call "reference calculations": namely calculating the energies and the forces with and without water in a hope to achieve stable forces calculation. This is what was suggested by you once on the APBS mailing list. In fact I used my pqr file (with and without 3rd line)  ATOM 1 I ION 1 3.000 0.000 0.000 1.000 1.000 ATOM 2 I ION 1 2.000 0.000 0.000 1.000 1.000 ATOM 2 I ION 1 2.000 2.000 0.000 1.000 1.000  for both mgauto and mgmanual (see the attached apbs.in file) APBS calculations. The attached config contains two calculations, with sdie=pdie=1 and pdie=1 and sdie=80. At the end I request calculateforce 1 (calculate total forces) and "print forces 2  1 end" in a hope to get the difference (the composition of the APBS infile is essentially taken from your examples) When two ions (the first two lines of the complex.pqr) are set for the calculation, the results are wonderful: Global net fixed charge force = (2.389761176097E05, 4.436307650468E05, 4.584334501816E05) kJ/mol/A If the third ion is added, I get Global net fixed charge force = (3.934905783681E+01, 1.583010217317E+02, 9.327207303476E05) kJ/mol/A Global net ionic boundary force = (0.000000000000E+00, 0.000000000000E+00, 0.000000000000E+00) kJ/mol/A Global net dielectric boundary force = (1.414546118941E+02, 2.060577689141E+02, 3.208009916572E06) kJ/mol/A Global net apolar boundary force = (2.626833777460E02, 3.726321592141E02, 2.098239783487E18) kJ/mol/A The resulting force is huge and does not seem to be eliminated by the subtraction. Adding the results from the coulomb binary would not help, since being a result of pair summation, the coulomb binary outputs the forces automatically respecting the 3rd newton law. The value of the forces I get do not depend strongly on the *len parameters. It seems like if I take other complexes from your examples and try to calculate total forces, I get large nonzero total forces. Does it mean that a) I am doing something stupid and wrong or b) the forces are not very reliable in this kind of calculations. Do I understand the APBS outputs right (namely I expect total forces printed out with the “print force ... end” statement?) I have tried also to introduce Debye screening to make sure my boundaries are far away enough to get "pure pair interaction". Did not help. Can you please suggest me something? I am playing with the program for quite awhile already and I did spend time reading the manuals. I am sure you are getting a number of questions like this. If there is something that helps me, I can compile our correspondence into a text for your site and give it back to you as a case study. I am sure others will be interested. To make my point maybe more clear: the example with two ions works well both in my settings and in your provided examples. Adding the third ion leads to crazy results (regarding the forces at least). There is another (almost offtopic) question: can you recommend a good general purpose nonlinear multigrid solver? Thank you very much for your time and good answer (and your software!), will be looking forward to hear from you again. Peter __________________________________________________________________ Introducing the New Netscape Internet Service. Only $9.95 a month  Sign up today at http://isp.netscape.com/register Netscape. Just the Net You Need. 
From: Nathan A. Baker <baker@ch...>  20040225 18:09:47

This is very strange. I just checked the force code and it is correct; i.e., the forces returned by APBS give the same results as the manual differentiation of the energy with respect to atom position. I'll continue to look into this problem and let you know what I find. Thanks for pointing this out.  Nathan Peter Fedichev <PeterOFedichev@...> (02252004 05:12:070500): >Dear Nathan, > >thank you very much for your mail. I have looked through it carefully and still would like to ask you for further clarification. > >I do understand that any numerical solver is going to have selfenergies and selfforces problems. This is why I am aiming at the forces first rather than for the energies. I also know that pointlike charges are fundamentally problematic for multigrid methods. > >With all this in mind I did try to do what you call "reference calculations": namely calculating the energies and the forces with and without water in a hope to achieve stable forces calculation. This is what was suggested by you once on the APBS mailing list. > >In fact I used my pqr file (with and without 3rd line) > >ATOM 1 I ION 1 3.000 0.000 0.000 1.000 1.000 >ATOM 2 I ION 1 2.000 0.000 0.000 1.000 1.000 >ATOM 2 I ION 1 2.000 2.000 0.000 1.000 1.000 > > >for both mgauto and mgmanual (see the attached apbs.in file) APBS calculations. The attached config contains two calculations, with sdie=pdie=1 and pdie=1 and sdie=80. At the end I request calculateforce 1 (calculate total forces) and "print forces 2  1 end" in a hope to get the difference (the composition of the APBS infile is essentially taken from your examples) > >When two ions (the first two lines of the complex.pqr) are set for the calculation, the results are wonderful: > >Global net fixed charge force = (2.389761176097E05, 4.436307650468E05, 4.584334501816E05) kJ/mol/A > >If the third ion is added, I get > Global net fixed charge force = (3.934905783681E+01, 1.583010217317E+02, 9.327207303476E05) kJ/mol/A > Global net ionic boundary force = (0.000000000000E+00, 0.000000000000E+00, 0.000000000000E+00) kJ/mol/A > Global net dielectric boundary force = (1.414546118941E+02, 2.060577689141E+02, 3.208009916572E06) kJ/mol/A > Global net apolar boundary force = (2.626833777460E02, 3.726321592141E02, 2.098239783487E18) kJ/mol/A > >The resulting force is huge and does not seem to be eliminated by the subtraction. Adding the results from the coulomb binary would not help, since being a result of pair summation, the coulomb binary outputs the forces automatically respecting the 3rd newton law. The value of the forces I get do not depend strongly on the *len parameters. It seems like if I take other complexes from your examples and try to calculate total forces, I get large nonzero total forces. > >Does it mean that a) I am doing something stupid and wrong or b) the forces are not very reliable in this kind of calculations. Do I understand the APBS outputs right (namely I expect total forces printed out with the ???print force ... end??? statement?) > >I have tried also to introduce Debye screening to make sure my boundaries are far away enough to get "pure pair interaction". Did not help. > >Can you please suggest me something? I am playing with the program for quite awhile already and I did spend time reading the manuals. I am sure you are getting a number of questions like this. If there is something that helps me, I can compile our correspondence into a text for your site and give it back to you as a case study. I am sure others will be interested. > >To make my point maybe more clear: the example with two ions works well both in my settings and in your provided examples. Adding the third ion leads to crazy results (regarding the forces at least). > >There is another (almost offtopic) question: can you recommend a good general purpose nonlinear multigrid solver? > >Thank you very much for your time and good answer (and your software!), >will be looking forward to hear from you again. > >Peter > >__________________________________________________________________ >Introducing the New Netscape Internet Service. >Only $9.95 a month  Sign up today at http://isp.netscape.com/register > >Netscape. Just the Net You Need. End of message from Peter Fedichev .  Nathan A. Baker, Assistant Professor Washington University in St. Louis School of Medicine Dept. of Biochemistry and Molecular Biophysics Center for Computational Biology 700 S. Euclid Ave., Campus Box 8036, St. Louis, MO 63110 Phone: (314) 3622040, Fax: (314) 3620234 URL: http://www.biochem.wustl.edu/~baker PGP key: http://cholla.wustl.edu/~baker/pubkey.asc 
From: Nathan A. Baker <baker@ch...>  20040225 18:58:35

Wow  this is very disturbing indeed. Here's what I've found so far: * The nonzero net force is independent of direction (I've tried it with the x and y axes flipped) * The nonzero component of the net force is *VERY* dependent on the grid spacing: SPACING NONZERO COMPONENT MAGNITUDE   0.156 1.598970527865E+01 0.104 3.719502847071E+01 0.078 6.079635610502E+00 0.062 1.273769475566E03 and grid position (at a fixed 0.078 spacing): XCENTER NONZERO COMPONENT MAGNITUDE   0.00 1.793940864605E+01 0.05 1.544555563973E+01 0.10 2.625380227759E+00 0.15 1.982065110814E+01 * The guilty component appears to be the dielectric boundary force. * For reasonable grid spacings, the net force seems to be largely independent of solute size. This does not appear to be a bug in the code; instead it is a consequence of the relatively sharp boundaries used in the PB definition. The upshot of all this is that one should use the finest grid possible when calculating PBbased forces. Peter Fedichev <PeterOFedichev@...> (02252004 05:12:070500): >Dear Nathan, > >thank you very much for your mail. I have looked through it carefully and still would like to ask you for further clarification. > >I do understand that any numerical solver is going to have selfenergies and selfforces problems. This is why I am aiming at the forces first rather than for the energies. I also know that pointlike charges are fundamentally problematic for multigrid methods. > >With all this in mind I did try to do what you call "reference calculations": namely calculating the energies and the forces with and without water in a hope to achieve stable forces calculation. This is what was suggested by you once on the APBS mailing list. > >In fact I used my pqr file (with and without 3rd line) > >ATOM 1 I ION 1 3.000 0.000 0.000 1.000 1.000 >ATOM 2 I ION 1 2.000 0.000 0.000 1.000 1.000 >ATOM 2 I ION 1 2.000 2.000 0.000 1.000 1.000 > > >for both mgauto and mgmanual (see the attached apbs.in file) APBS calculations. The attached config contains two calculations, with sdie=pdie=1 and pdie=1 and sdie=80. At the end I request calculateforce 1 (calculate total forces) and "print forces 2  1 end" in a hope to get the difference (the composition of the APBS infile is essentially taken from your examples) > >When two ions (the first two lines of the complex.pqr) are set for the calculation, the results are wonderful: > >Global net fixed charge force = (2.389761176097E05, 4.436307650468E05, 4.584334501816E05) kJ/mol/A > >If the third ion is added, I get > Global net fixed charge force = (3.934905783681E+01, 1.583010217317E+02, 9.327207303476E05) kJ/mol/A > Global net ionic boundary force = (0.000000000000E+00, 0.000000000000E+00, 0.000000000000E+00) kJ/mol/A > Global net dielectric boundary force = (1.414546118941E+02, 2.060577689141E+02, 3.208009916572E06) kJ/mol/A > Global net apolar boundary force = (2.626833777460E02, 3.726321592141E02, 2.098239783487E18) kJ/mol/A > >The resulting force is huge and does not seem to be eliminated by the subtraction. Adding the results from the coulomb binary would not help, since being a result of pair summation, the coulomb binary outputs the forces automatically respecting the 3rd newton law. The value of the forces I get do not depend strongly on the *len parameters. It seems like if I take other complexes from your examples and try to calculate total forces, I get large nonzero total forces. > >Does it mean that a) I am doing something stupid and wrong or b) the forces are not very reliable in this kind of calculations. Do I understand the APBS outputs right (namely I expect total forces printed out with the ???print force ... end??? statement?) > >I have tried also to introduce Debye screening to make sure my boundaries are far away enough to get "pure pair interaction". Did not help. > >Can you please suggest me something? I am playing with the program for quite awhile already and I did spend time reading the manuals. I am sure you are getting a number of questions like this. If there is something that helps me, I can compile our correspondence into a text for your site and give it back to you as a case study. I am sure others will be interested. > >To make my point maybe more clear: the example with two ions works well both in my settings and in your provided examples. Adding the third ion leads to crazy results (regarding the forces at least). > >There is another (almost offtopic) question: can you recommend a good general purpose nonlinear multigrid solver? > >Thank you very much for your time and good answer (and your software!), >will be looking forward to hear from you again. > >Peter > >__________________________________________________________________ >Introducing the New Netscape Internet Service. >Only $9.95 a month  Sign up today at http://isp.netscape.com/register > >Netscape. Just the Net You Need. End of message from Peter Fedichev .  Nathan A. Baker, Assistant Professor Washington University in St. Louis School of Medicine Dept. of Biochemistry and Molecular Biophysics Center for Computational Biology 700 S. Euclid Ave., Campus Box 8036, St. Louis, MO 63110 Phone: (314) 3622040, Fax: (314) 3620234 URL: http://www.biochem.wustl.edu/~baker PGP key: http://cholla.wustl.edu/~baker/pubkey.asc 