Hi there,

I’m writing back with an update - we tried a few modifications to the procedure I described to derive the sialic acid charges.  The “original” procedure is described earlier in the thread, and each test involves changing only one thing in the whole procedure.  Here are some RMS errors, comparing the charges we derived vs. the charges in GLYCAM.

- RMSE vs. GLYCAM charges -
Our “original” procedure: 0.0572
Ring dihedrals constrained (in addition to exocyclic) in the geometry optimization step: 0.0580
Adding diffuse functions (6-31+G*) in the geometry optimization and ESP calculation steps: 0.0807
Using ROH as the aglycon (instead of OME): 0.0550

From these tests, I think that adding diffuse functions actually increases the difference vs. the GLYCAM charges, and the other changes only make a small difference.  I’m still not sure what is preventing us from reproducing the GLYCAM charges more closely, but perhaps this level of error is acceptable for proceeding with our research.  Please let me know if you have additional ideas.

Sincerely,

- Lee-Ping

On Jul 26, 2018, at 6:11 PM, Lee-Ping Wang <[log in to unmask]> wrote:

Hello Lachele,

Thank you very much for your feedback.  I was traveling for a few days, sorry for the slow response.  The large distortions in the ring (after constrained energy minimization at HF/6-31G*) are quite common for this system.  I expect they may become even more common if we add more chemical substitutions.  

I will try two approaches and get back to you with results: (1) Constraining the ring dihedral angles in addition to the exocyclic ones when doing energy minimization; and (2) including diffuse functions in the basis set when minimizing and computing the ESP.

Sincerely,

- Lee-Ping

On Tue, Jul 24, 2018 at 10:33 PM Lachele Foley <[log in to unmask]> wrote:
In the course of an MD simulation, it is common to encounter structures that are not at a local minimum.  If a structure is flexible - that is, if its energy landscape is 'flatter' - then it might stray quite far from local minima.  I am not familiar with the finer details of sialylate's ring-puckering prefeences, but its electronic structure causes it to have other interesting conformational proclivities, notably increased flexibility in the psi dihedral (the exo-anomeric effect).

In the image you show, I suspect that the main drivers for the conformational change are attraction between the H5N and the carboxylate and repulsion between the regions of high ESP in the O5N and O4.  In this case, it might be more realistic to find a way to sample conformations more similar to the un-minimized structure.  This is especially so if structures such as these account for a reasonable fraction of the structures from your simulations.  If these changes are relatively rare, it might not make a lot of difference.

Diffuse functions are needed for proper charge calculations in charged species.  I don't know enough to say how strong their effect is on geometries in these systems.  You should take a look at the papers by Matt Tessier (for examples, see link below).  He did a lot of the work with NeuNAc.  Others in the group might also have opinions.


On Tue, Jul 24, 2018 at 12:48 AM Lee-Ping Wang <[log in to unmask]> wrote:
Hello Xiaocong,

The use of different basis sets is one possible source of the difference.  I based the choice of using 6-31G* on the following:

From the 2007 GLYCAM paper, there is the quote “For each of these snapshots partial charges were calculated by fitting to the HF/6-31G* MEP. Prior to the charge calculations, each structure was optimized at the HF/6-31G* level, with the rotatable exocyclic bonds constrained to their MD conformations.” In addition, the RED server also uses this basis set to derive GLYCAM charges:

# RESP-C2:  HF/6-31G(d)//HF/6-31G(d)                              CHELPG algo.    1 RESP fit(*1)  qwt=.0100
#   Used in the Glycam FF

On the other hand, the 2007 paper also says "The HF/6-31G* level of theory was employed for neutral fragments, whereas for anionic molecules diffuse functions were added.”  The basis set you mentioned, 6-311+G**, differs from 6-31G* in several ways – not only does it add diffuse functions, but also a valence shell on the heavy atoms and polarization functions on the hydrogens.  Is there a better reference for the charge derivation procedure that you could point me to?

Thanks a lot,

- Lee-Ping

On Jul 23, 2018, at 8:33 PM, Xiaocong Wang <[log in to unmask]> wrote:

Hi Dr. Wang,
We usually have HF/6-311+G** for sugars with charges (both geometry optimization and ESP calculations).  I am not sure if this would be the problem. 
I have the same thoughts for where the problems might also be.  I usually don’t see ring distortion like that.  But, I did not derive charges for sialic acid.  I will ask around for people in our lab who have worked on this sugar to see if the ring distortion in it is common.
Thank you!
Best regards,
Xiao


Xiaocong Wang
Complex Carbohydrate Research Center
The University of Georgia
315 Riverbend Road,
Athens, GA, 30602
Tel: (706) 254-7958

From: Users of GLYCAM & GLYCAM-Web <[log in to unmask]> on behalf of Lee-Ping Wang <[log in to unmask]>
Sent: Monday, July 23, 2018 11:24:15 PM
To: [log in to unmask]
Subject: Re: Comparing self-derived GLYCAM charge parameters
 
Hello Xiaocong,

Thanks for your reply!

1) We used HF/6-31G* for constrained energy minimization and also for the ESP calculations.  

2) We set the charges on certain hydrogens in 0SA to zero, based on whether they were bonded to carbon. The atom names are: H3E, H3A, H4, H5, H6, H7, H8, H9S, H9R, H3M, H2M, H1M. I think these are the same hydrogens that got zero charges in GLYCAM. For OME we froze all of the charges and did not optimize them at all, which meant hydrogens H1, H2, H3 had charges of zero in the calculation.

At the moment, the only “problem” I can imagine is that the constrained energy minimizations might be causing some unexpected large distortions of the sugar ring, because the ring dihedral angles are not being constrained (e.g. see the picture highlighting C3-C4-C5-C6 from the original email).  I am attaching a table of our ensemble-averaged charges and the standard deviations, and the standard deviations of carbon charges are rather large (σ = 0.05 ~ 0.13).  

Have you ever seen major distortions of the 6-membered ring in the constrained minimizations?  Do you think that constraining the ring dihedral angles is a good idea?

Thanks,

- Lee-Ping



On Jul 23, 2018, at 6:57 PM, Xiaocong Wang <[log in to unmask]> wrote:

Dear Dr. Wang,
Your procedure is the same as what we did for deriving ensemble-averaged partial atomic charge sets for carbohydrate molecules.  It is a very thorough work for reproducing charge sets for 0SA.  I have two questions about calculation detail:  1. What level of theories did you use for geometry optimization and ESP calculation? 2. In step 6, by "aliphatic hydrogens", you mean all aliphatic hydrogen atoms in 0SA+OME, not just those in OME?
The purpose of ensemble-averaged charge sets is to have better representations of electrostatic properties for exocyclic groups.  From my experience, charge values for ring carbon atoms are usually stable.  I am not sure what caused this variation, but I am happy to take a closer look if I can get more info.
Best regards,
Xiao

Xiaocong Wang
Complex Carbohydrate Research Center
The University of Georgia
315 Riverbend Road,
Athens, GA, 30602
Tel: (706) 254-7958

From: Users of GLYCAM & GLYCAM-Web <[log in to unmask]> on behalf of Lee-Ping Wang <[log in to unmask]>
Sent: Monday, July 23, 2018 6:29:17 PM
To: [log in to unmask]
Subject: Comparing self-derived GLYCAM charge parameters
 
Hello there,

We’re carrying out some studies on chemically modified sialic acids - specifically N-acetylation at the C4, C7, C8, or C9 positions.  Because there are many potential acetylation sites, we became interested in questions such as: (1) how much do the charges depend on the N-acetylation site, and (2) how much does the N-acetylation affect the charges on the sialic acid itself.  

For the initial step in our studies, we tried to reproduce your charge derivation for unmodified sialic acid (residue 0SA).  We obtained results that were quite close to the original GLYCAM charges, but there are still some differences and I’m not sure how significant they are.  I’d like to know if we’re following your procedure correctly, or if there’s anything we should be doing differently to improve the agreement.

Here is the procedure we used:

1) Using tleap, we created the oligosaccharide sequence OME 0SA (using the GLYCAM_06j-1 force field).
2) The system was neutralized by adding one Na+ and solvated in a box of TIP3P water.
3) We ran 200 ns of unbiased MD simulation and saved 100 structures of the oligosaccharide equally spaced at 2 ns intervals.
4) We energy-minimized each structure at the HF/6-31G* level of theory with all exocyclic dihedral angles constrained to their values from the MD structure.  (Note: these calculations sometimes cause the ring dihedral angles to change significantly.  In one calculation, the C3-C4-C5-C6 dihedral angle increased from -66 to +13 degrees; see picture.)
5) We generated grids and calculated ESP values for each minimized structure using Gaussian and the options "Pop=chelpg  IOp(6/33=2)”.
6) We ran RESP on each individual structure using the resp program from AmberTools.  
Charges on all atoms were restrained (including hydrogen) with a weight of 0.01.
Charges on OME were held constant and charges on aliphatic hydrogens were set to zero.  
We preferred to use the old resp program rather than RED server because the former gave us a better understanding of the process.
7) The set of 100 RESP-optimized charges were averaged to obtain the final result.
8) The final charges are plotted vs. the GLYCAM charges for comparison (see picture).

The largest charge difference was on C2: 0.237 (GLYCAM) vs. 0.069 (our charges), a difference of 0.168.  
For polar atoms (|q| > 0.5) the largest difference was on C1: 0.789 (GLYCAM), vs. 0.874 (our charges), a difference of 0.087.
Generally, the most significant differences were on the carbon atoms (+/- 0.04–0.10 for C3–C9).

9) We also repeated steps 3-8 with a 500 ns simulation, taking 100 snapshots equally spaced at 5 ns intervals.  The results are very close to our 200 ns results, with the greatest charge difference being 0.013.  With a 1000 ns sampling simulation, the largest charge difference vs. 200 ns is 0.023.  Thus, I do not think our results depend much on the trajectory length.

I’m very interested in your thoughts on these results.  Thanks again!

Sincerely,

Lee-Ping Wang
Assistant Professor
Department of Chemistry, UC Davis

<Geometries.png>

<Charges.png>



--
:-) Lachele
Lachele Foley
CCRC/UGA
Athens, GA USA