Previous Thread
Next Thread
Print Thread
Page 1 of 2 1 2
Dihedral with k=0
#33513 03/03/14 03:11 PM
Joined: Oct 2010
Posts: 63
C
chemist Offline OP
Forum Member
OP Offline
Forum Member
C
Joined: Oct 2010
Posts: 63
Hello, I'm trying to parametrise the DOTA chelator: 1,4,7,10-tetraazacyclododecane-1,4,7,10-tetraacetic acid with the formula (CH2CH2NCH2CO2H)4. This I have attached to a standard amino acid (I've chosen Lysine).

The parameters I need to optimise all involve the DOTA tertiary amine junction to the Lysine branch. Not many cgenff parameters exist for tertiary amines.

This is a big molecule (73 atoms) so I've split the DOTA part (not Lys) into 2 fragments.

One particular dihedral is causing problem to fit. The penalty is below 50, but it must be validated:

NG2S1 CG2O1 CG321 NG301 0.4000 1 0.00 ! Kw_mol , from NG2S1 CG2O1 CG324 NG3P3, penalty= 37.9

OG2D1 CG2O1 CG321 NG301 0.0000 1 0.00 ! Kw_mol , from OG2D1 CG2O1 CG324 NG3P3, penalty= 37.9

CG2O1 CG321 NG301 CG321 2.5000 1 180.00 ! Kw_mol , from NG311 CG321 NG311 CG2R61, penalty= 146.5
CG2O1 CG321 NG301 CG321 1.5000 2 0.00 ! Kw_mol , from NG311 CG321 NG311 CG2R61, penalty= 146.5
CG2O1 CG321 NG301 CG321 0.5000 3 0.00 ! Kw_mol , from NG311 CG321 NG311 CG2R61, penalty= 146.5


The parameter is derived from OG2D1 CG2O1 CG324 NG3P3 which also has n=1, k=0.000. These non-contributing dihedrals are next to a rigid peptide bond, so I'm guessing we do not want them contributing to the potential energy function as they would disrupt planarity?

Do I fit a bond like this? Do I use the Paramchem guess?

Here I tried fitting all 3 dihedrals, and then I will use the Paramchem values for the first two, and the fitting values for the n=1,n=2, n=3 dihedral: www.peecee.dk/upload/view/433133



Also: Kenno, could I ask for the parameters you have for PAP-1 coumarin?

Thanks

Last edited by chemist; 03/03/14 10:18 PM.
Re: Dihedral with k=0
chemist #33529 03/07/14 08:15 PM
Joined: Dec 2005
Posts: 1,535
Forum Member
Offline
Forum Member
Joined: Dec 2005
Posts: 1,535
First a general note: I understand you chose the neutral amine because it's chelating a metal. If so, you'll need to worry about representing the metal chelation correctly, which can be really tricky depending on the metal, its oxidation state and its coordination. Here's a related discussion. if not, you should protonate the nitrogen.

Originally Posted By: chemist
The parameter is derived from OG2D1 CG2O1 CG324 NG3P3 which also has n=1, k=0.000. These non-contributing dihedrals are next to a rigid peptide bond, so I'm guessing we do not want them contributing to the potential energy function as they would disrupt planarity?
It's complicated. The rotation around that bond is controlled by 4 parameters: OCCN, OCCC(beta), NCCN and NCCC(beta). The final profile is the sum of the 1-4 and longer range non-bonded interactions + all the multiplicities of all the dihedrals, each of which is either added or subtracted depending on the multiplicity and geometrical considerations. Also related. We generally limit our choice of multiplicities to the ones that are geometrically and chemically sensible. Also, we want our parameters to be transferable, which leads us to constrain ourselves even further, eg. by setting some types of dihedral parameters always to a predefined value and only fitting the remaining ones, or sometimes by requiring certain fitted parameters to have the same values. Finally, when dihedral contributions with high amplitudes partially or completely cancel out, this leads to a residual angular force which is generally undesirable. All these considerations are taken into account when fitting to the QM. Then, in the next step, bulk phase simulations are done on bigger relevant model systems, and the dihedrals are tweaked to match experimental (often NMR) data, again taking into account all the above considerations. So if you look at the end product, and ask: "why does this parameter have this value", the question is not always easy to answer. You are right that setting the OG2D1 CG2O1 CG324 NG3P3 to a high value and counteracting this with another parameter would disrupt the planarity, and that likely played a role in the decision to set it to 0, but a low value should not be disastrous.

Looking forward, what we would do in the most generic case is scan and fit the N-C-C-N because that's the backbone of your molecule, as opposed to the O, which is a dead end. Prior to that, N-C-C-O is fitted to the C=O wagging motions. However, this generic reasoning is not valid here for 2 reasons. First, we know that the C=O is part of a rigid amide group and that the bond currently under consideration is a free rotor, so it only makes sense to fit the dihedrals on the amide side to the wagging motions, not the dihedral you're talking about. Second, wagging motions of carbonyl groups that are not part of rigid planar rings cannot be fit by angles and dihedrals alone, so an improper is needed, which becomes the primary agent for getting the wagging motion right.

Are you still there? Yeah, force field parameterization is a lot more complicated than meets the eye.

Edit: oh yeah, you probably want a recommendation. I think keeping this one at 0 is not a bad idea. Especially if you're using automatic fitting; it will make things easier on the poor stupid fitting algorithm. (No offense to its creators; I'm calling all fitting algorithms stupid because they don't have common sense or chemical intuition. It's really really hard to put that stuff into a box.)

Also, you attached an energy plot, but such plots don't say much if you don't also show the parameters. If you have a bad plot and the parameters are weird, then it is likely something has gone wrong with the fitting. If you get a good plot and the parameters are weird, that's still not good; they may be overfitted and blow up in your face when transfered. Conversely, if you get a mediocre plot but the parameters are sensible, it might just be that the force field is too approximate to capture physical reality more accurately than that.

Originally Posted By: chemist
Also: Kenno, could I ask for the parameters you have for PAP-1 coumarin?
I finally posted them.

Last edited by Kenno; 03/07/14 10:48 PM.
Re: Dihedral with k=0
chemist #33616 03/25/14 03:07 PM
Joined: Oct 2010
Posts: 63
C
chemist Offline OP
Forum Member
OP Offline
Forum Member
C
Joined: Oct 2010
Posts: 63
Hello Kenno,

Thanks for this very detailed and clear answer.

I still have some questions about the parametrisation of the DOTA chelator in the presence of Gd(III) ion:

Quote:


If you really must, you can put some dummy parameters on the metal ion and use constraints to get the coordination geometry right, but even then, it will be tricky to get the electrostatics right (and you'd have to be careful to avoid destabilizing the MD simulation).




1. I wish to parametrise DOTA with the metal.

Does that mean I:
I) Carry out QM calculations (water interaction, Hessian, and dihedral scans) on the full system.
II) Perform the MM part with a dummy parameter for Gd?


There are CHARMM parameters for the LJ parameters of Gd(III). I can use those for the dummy atom.


2. If I use a dummy for the MM fitting part:

I set the charge of the dummy to a value which will balance out to an total integer charge. I have RESP values of the partial charges, as well as GD-DOTA bonded parameters published by Henriques et al., Int. Journal of Quantum Chemstry, 1999.

3. For the CHARMM parametrisation:
You advocate a splitting approach, but when there's a metal bound, how do I split the chelator? I can't see how this would be done.


4. For the actual MD simulation:

I have seen in the literature that a harmonic restraint approach is used on the dummy bond, angle and dihedral using e.g. the NAMD Colvars feature.

In my case, I can access 2 conformations according to NMR. Do I fit / restrain relative to the minimum energy conformer?






--------


Alternatively:

You write I should optimise the chelator with protonated N. Fine, but at the end of the procedure, I will have a topology file for protonated DOTA. How could I modify this to account for the metal?
I can think of 2 approaches:
I) Create a second topology file without N protonated. The charge imbalance is accounted for by the presence of a dummy metal parameter with charge balancing the 4 protons?


II) Parametrise the whole DOTA chelator partial charges to add up to zero when excluding the protons.
Then define a DOTA topology file
And define a GD topology file containing GD with integer +3 charge. Use bond patches to link DOTA residue with GD residue.
This ignores polarisation of charges.






Are you aware of any resources for metal parametrisation in CHARMM?



Thanks

Re: Dihedral with k=0
chemist #33620 03/25/14 04:41 PM
Joined: Dec 2005
Posts: 1,535
Forum Member
Offline
Forum Member
Joined: Dec 2005
Posts: 1,535
Originally Posted By: chemist
II) Perform the MM part with a dummy parameter for Gd?
No, I talked about the dummy parameters assuming you do not want to parameterize them metal. Which I thought was a safe assumption; you're really getting yourself into a hairy situation if you want to do that.

Originally Posted By: chemist
There are CHARMM parameters for the LJ parameters of Gd(III). I can use those for the dummy atom.
Not in the official release, there aren't. I know there are a lot of unofficial parameter sets for various metals out there, but their quality is very spotty, so I'd need to have a look at what kind of parameterization effort is behind it.


Originally Posted By: chemist
I set the charge of the dummy to a value which will balance out to an total integer charge. I have RESP values of the partial charges, as well as GD-DOTA bonded parameters published by Henriques et al., Int. Journal of Quantum Chemstry, 1999.
If you go for the pragmatic route and don't parameterize the metal (restraining its coordination instead), you can do pretty much anything. If you do want to parameterize the metal, then the charges you got are useless. Among our Frequently Asked Questions, this is a particularly frequent one.

Originally Posted By: chemist
You advocate a splitting approach, but when there's a metal bound, how do I split the chelator? I can't see how this would be done.
I'd need to have some idea of the 3D structure of the complex in order to answer that question (and to give more definite advice on whether there's any point at all in doing the full parameterization, or you should just go for the pragmatic route).

Originally Posted By: chemist
In my case, I can access 2 conformations according to NMR. Do I fit / restrain relative to the minimum energy conformer?
Same as above, but if they are both relevant, and differ to the extent they can't readily interconvert, then the best thing to do would probably be two independent simulations for the two species.

Originally Posted By: chemist
You write I should optimise the chelator with protonated N.
No, I wrote you should do that if not [chelating a metal]. Sorry if I formulated that paragraph ambiguously.

Originally Posted By: chemist
Are you aware of any resources for metal parametrisation in CHARMM?
For free metals in solution, Benoît Roux et al.'s work is authoritative. For metals that are bound in a complex, well, there have been various efforts with various (usually very high) degrees of approximation, but there's very little consistency, and nothing stands out as a "gold standard". Nothing except QM/MM, that is.

Re: Dihedral with k=0
chemist #33621 03/25/14 04:46 PM
Joined: Oct 2010
Posts: 63
C
chemist Offline OP
Forum Member
OP Offline
Forum Member
C
Joined: Oct 2010
Posts: 63
Hello Kenno,

My top priority is the ligand conformation, so let's say I put the metal parameters as lowest priority, I would use a dummy.


Could you explain how exactly you would parametrise that system? If I don't parametrise the metal, do I put a dummy during the parametrisation?

Can I give the dummy a formal charge?

For the quantum calculations, I assign the dummy a type X.
For the MM fitting, what charge sum will I aim for the chelator? The opposite value of the formal charge of the dummy metal?

Thanks

Re: Dihedral with k=0
chemist #33623 03/25/14 04:52 PM
Joined: Oct 2010
Posts: 63
C
chemist Offline OP
Forum Member
OP Offline
Forum Member
C
Joined: Oct 2010
Posts: 63
I try to clarify: What's the chemistry of the system I would parametrise?

since it is DOTA-GD, if I don't parametrise the metal, then would I be parametrising DOTA-DZ, where DZ is my dummy.



http://upload.wikimedia.org/wikipedia/commons/thumb/f/f8/Gadoteric_acid.png/220px-Gadoteric_acid.png

As you see, splitting this is extremely tricky if a metal or a dummy is bound.

Re: Dihedral with k=0
chemist #33624 03/25/14 04:54 PM
Joined: Dec 2005
Posts: 1,535
Forum Member
Offline
Forum Member
Joined: Dec 2005
Posts: 1,535
  • If your ligand is tightly bound to the metal, you can't really consider it independently from the metal, especially in the QM calculations.
  • I wrote "dummy parameters" not "dummy atom". Using a dummy atom for your metal in the QM will have the strangest results.
  • I realize what I wrote so far may seem contradictory; I really feel I cannot give sensible advice without more information.

Re: Dihedral with k=0
chemist #33626 03/25/14 05:02 PM
Joined: Dec 2005
Posts: 1,535
Forum Member
Offline
Forum Member
Joined: Dec 2005
Posts: 1,535
OK, our posts crossed each other. Going by the 2D structure you posted, you'll definitely have to follow the pragmatic route, but some details still need to be filled in. I have to sleep over it for a night, and at any rate, my self-imposed CHARMM forum time quotum is exceeded for today. Remind me if I don't reply in the next 36h or so. In the meanwhile, it would be really helpful if you could give me some 3D information on the 2 conformations.

Re: Dihedral with k=0
chemist #33648 03/26/14 12:07 PM
Joined: Oct 2010
Posts: 63
C
chemist Offline OP
Forum Member
OP Offline
Forum Member
C
Joined: Oct 2010
Posts: 63
Hello Kenno,

Firstly, I must clarify I'm working on a branched version of DOTA, with a lysine branch.

For me, this parametrisation is extremely urgent, so I prefer to work on the fastest method available. We put as lowest priority the accurate representation of the metal, which is poorly treated in Class I ff anyway.


DOTA-GD comes in 2 geometries: Square antiprismatic and Inverse square antiprismatic. These differ in the orientation (clockwise or counterclockwise) of the carboxylate branches.

so from what I understand you advise:

I) Run QM calculations on DOTA-GD, using the relevant basis sets to converge the metal correctly.

# MP2(FULL)/GEN 6D 7F Opt=(Redundant) SCF=Tight pseudo=read Geom=PrintInputOrient



II) Perform MM fit on the chelator. Use dummy parameters on metal for fitting?

III) Run MD with restrained dummy parameters on metal.




I upload the square antiprismatic (Lys-branched) DOTA-GD chelator.

http://peecee.dk/upload/view/434778
http://peecee.dk/upload/view/434779
http://peecee.dk/upload/view/434780

Attached Files
DOTA.png (85.11 KB, 388 downloads)
Re: Dihedral with k=0
chemist #33649 03/26/14 12:10 PM
Joined: Oct 2010
Posts: 63
C
chemist Offline OP
Forum Member
OP Offline
Forum Member
C
Joined: Oct 2010
Posts: 63
Or alternative approach:


I) Parametrising chelator in the absence of metal (not in PDB/PSF)

II) Adding restrained dummy parameters to mimic the metal during the MD

Re: Dihedral with k=0
chemist #33655 03/26/14 03:47 PM
Joined: Oct 2010
Posts: 63
C
chemist Offline OP
Forum Member
OP Offline
Forum Member
C
Joined: Oct 2010
Posts: 63
So my question is: What would be the simplest way of accounting for the presence of the metal in calculating the partial charges of the system?


When I submitted the chelator in Paramchem, I am getting charges that do not account for the presence of a metal.

Would submitting the protonated form of DOTA (all 4 N protonated) represent a reasonable approximation of the polarised charges by the metal?

Last edited by chemist; 03/26/14 03:47 PM.
Re: Dihedral with k=0
chemist #33670 03/27/14 03:42 PM
Joined: Oct 2010
Posts: 63
C
chemist Offline OP
Forum Member
OP Offline
Forum Member
C
Joined: Oct 2010
Posts: 63
To update you:

Going for the pragmatic route:

I took the CHARMM parameters I generated for the lysine-branched chelator (charges, bonds, angles and dihedrals), and inserted the metal assigning it a dummy charge of +3.
I restrained only the metal-chelator bonds and angles to values that were generated to reproduce experimental NMR values in the paper:

Computational approaches to the study of some lanthanide (III)–polyazamacrocyclic chelates for magnetic resonance imaging

Elsa S. Henriques, Margarida Bastos, Carlos F. G. C. Geraldes, Maria João Ramos

International Journal of Quantum Chemistry (Impact Factor: 1.31). 03/1999; 73(2):237 - 248. DOI:10.1002/(SICI)1097-461X(1999)73:2<237::AID-QUA17>3.0.CO;2-W





Now, in that paper they minimise both conformers with the parameters and demonstrate a 3 kcal/mol energy difference, consistent with experiments. This, to them, is a way of deciding on the validity [sic] of their parametrisation. I know we're far away from the CHARMM philosophy here.

Re: Dihedral with k=0
chemist #33692 04/01/14 04:35 PM
Joined: Dec 2005
Posts: 1,535
Forum Member
Offline
Forum Member
Joined: Dec 2005
Posts: 1,535
I finally found some time to look over the structures and download and read the paper:
  • One of the reasons why I was so adamant on 3D structures was that I wanted to know if there was any chance that part of the metal would be solvent-exposed. Based on the current information, the answer to that question seem to be "reversible binding of the water to the metal can and does occur". This largely invalidates the proposal of using dummy L-J parameters and charges.
  • Although the methodology in the paper is different from our standard recommendations, we already established that the standard approach doesn't apply to this system for several reasons. Keeping this in mind, Henriques et al.'s approach seems to be quite solid, and I would recommend you just transfer as much as you can from their paper. The only thing I would be careful of is the "capping water"; since it has such exceptional charges, I would introduce a weak restraint (using CHARMM's NOE functionality) that is zero in the range of normal interaction distances and comes into effect when the water starts dissociating. The idea is to prevent that abnormal water from drifting off into the solvent and wrecking havoc.

For the reference of other people who come across this thread, I would like to repeat that I almost always recommend against using 3rd party parameter sets from the literature, especially if they don't follow standard procedures. It's just that for this very exceptional case, the standard procedure doesn't apply, and the literature approach is of a rarely seen quality. In fact, for this system, I cannot come up with a significantly better alternative approach than what Henriques et al. did (bar a 6-month parametrization project).

Edit: oh, for the record, the Henriques et al. charges do show some clear signs of overfitting and the L-J parameters would qualify as "poorly validated", but for a system like this, both flaws are very hard to avoid, and since it's conformationally quite rigid and we're not transfering the charges to a significantly different molecule, I think we can get away with it for this once.

Last edited by Kenno; 04/01/14 04:42 PM.
Re: Dihedral with k=0
chemist #33693 04/01/14 07:02 PM
Joined: Oct 2010
Posts: 63
C
chemist Offline OP
Forum Member
OP Offline
Forum Member
C
Joined: Oct 2010
Posts: 63
Hello Kenno,

ok, I had initially tried to mix their parameters with mine, but I was not confident it would work.

I found that my CHARMM parameters with restraints on the metal, a +3 charge on the metal, and the LJ parameters from Henriques, were able to reproduce the crystal structure of DOTA to a good extent (no data shown, I'm still working on processing all the data). The backbone RMSD of the DOTA in solution is 0.5 A rel. to the crystal, and I was otherwise going to consider it sufficient to start simulating.

The big problem with the Henriques paper is that they do not give parameters except for the Metal - Ligand interactions (so, bonds and angles to the metal). The partial charges are not at all clear how to use (no sign given).

I have found more sucess with this approach than with the scrambled mix I got by taking those Henriques et al. parameters and mixing them with mine...

Last edited by chemist; 04/01/14 07:10 PM.
Re: Dihedral with k=0
chemist #33694 04/01/14 07:30 PM
Joined: Dec 2005
Posts: 1,535
Forum Member
Offline
Forum Member
Joined: Dec 2005
Posts: 1,535
I don't think the +3 charge on the metal is very realistic. Highly charged entities transfer a lot of their charge to their immediate surroundings in bulk phase, and this has a relevant effect on how these immediate surrounding interact with other stuff. Also, what do you mean with "no sign given"? I'm looking at the RESP columns in table IV and don't see anything that would preclude usage.

As for the bonded parameters, they probably transfered a lot from the CHARMM22 protein FF. One thing you could do is e-mail the corresponding author asking for the force field files (just google her; she's still leading a research group). That would leave a lot of guesswork out of the process, and explaining out of your future papers.

Re: Dihedral with k=0
chemist #33695 04/01/14 09:23 PM
Joined: Oct 2010
Posts: 63
C
chemist Offline OP
Forum Member
OP Offline
Forum Member
C
Joined: Oct 2010
Posts: 63
Hello,

Well, I tried adding up the RESP column, to me it is not clear how I would achieve an overall integer value with their charges (I can't figure out which ones are + and which ones are -).

I already contacted her a while ago, this paper is about 15 years old, I doubt they still have the data around.


If the CHARMM parameters for the chelator with a restraint on the metal can yield an ensemble of structure in the bulk that reproduces the crystal structure, is this not better than mixing parameters? I have no guarantee the author still has the values they used with RESP (which is from the AMBER philosophy?).

Anyhow, thanks for your help, much appreciated!


Otherwise, I will have to give the CHARMM values another try.

Re: Dihedral with k=0
chemist #33703 04/02/14 06:39 PM
Joined: Dec 2005
Posts: 1,535
Forum Member
Offline
Forum Member
Joined: Dec 2005
Posts: 1,535
I think there might be something wrong with the pdf viewer you're using or perhaps with the fonts you have installed on your system; as you can see in the attached crop from a screenshot, I'm seeing positive and negative values.

If you cannot reach the group leader, maybe try the first author; group leaders are often very busy. The fact that it's 15 years ago is not necessarily a problem; I still have copies of my very first molecular modeling work, dated 1999-2000.

Originally Posted By: chemist
If the CHARMM parameters for the chelator with a restraint on the metal can yield an ensemble of structure in the bulk that reproduces the crystal structure, is this not better than mixing parameters?
I'm more concerned about how it interacts with water, which they seem to have done some effort to get right.

Originally Posted By: chemist
RESP (which is from the AMBER philosophy?).
The big problem with this molecule is the charges. Both RESP and the CHARMM-style water interactions have issues with buried atoms. In CHARMM, we work around these issues by splitting up the system, and then transfering charges when re-assembling, but here, you cannot do that because the +3 ion in the middle that is not completely shielded from solvent changes everything. If it would be possible to do a large number of water interactions, then it would in theory (though often not in practice) be possible to resolve the effect of the buried charges. However, the CHARMM methodology calls for constraining oneself to relevant water interactions, of which there are not enough. In this respect, the many ESP samples used in RESP have a definite advantage. Yes, the charges will not be balanced against the other parameters. However, AMBER also uses TIP3P water and the Lorentz-Berthelot combining rules, so they likely won't be totally outrageously catastrophical (compared to more exotic nonbonded models like MMFF94). Usually, I vehemently oppose using RESP charges (or the Merz-Kollman charges we used to use as an initial guess before the CGenFF program was ready) because there's always a way to get charges that are more likely to be CHARMM-compatible with minimal to moderate effort. However, in this particular case case, I don't see how I would do that, so the (likely poorly balanced) charges will have to do in a pinch. It can't possibly be worse than setting the metal to +3 and using transfered chrges for the atoms that bind to it (and are also solvent-exposed).

Attached Files evince_crop.png
Last edited by Kenno; 04/02/14 06:49 PM. Reason: added attachment, reference to Merz-Kollman, and last sentence.
Re: Dihedral with k=0
chemist #33709 04/02/14 08:09 PM
Joined: Oct 2010
Posts: 63
C
chemist Offline OP
Forum Member
OP Offline
Forum Member
C
Joined: Oct 2010
Posts: 63
Thanks, I will let you know when I have an update.


Yes, there was a problem with the PDF viewer indeed. I reopened it in Adobe outside the browser and the signs appear. The printed manuscript was also without the charges.


Chemist

Re: Dihedral with k=0
chemist #33710 04/02/14 08:15 PM
Joined: Dec 2005
Posts: 1,535
Forum Member
Offline
Forum Member
Joined: Dec 2005
Posts: 1,535
Is this by any chance Firefox' built-in pdf viewer? Because I don't like that one. Apart from occasional display issues, its performance is horrible. The only good thing about it is that it gets rid of the gaping security hole that is the Adobe PDF plugin.

Re: Dihedral with k=0
chemist #33712 04/02/14 09:13 PM
Joined: Oct 2010
Posts: 63
C
chemist Offline OP
Forum Member
OP Offline
Forum Member
C
Joined: Oct 2010
Posts: 63
I opened it in both the pdf viewer of Chrome and Firefox.

I could only see them when loading Adobe externally. It's weird.

Thanks.

C

Page 1 of 2 1 2

Moderated by  alex, lennart, rmv 

Link Copied to Clipboard
Powered by UBB.threads™ PHP Forum Software 7.7.4
(Release build 20200307)
Responsive Width:

PHP: 5.6.33-0+deb8u1 Page Time: 0.020s Queries: 56 (0.011s) Memory: 1.1026 MB (Peak: 1.3242 MB) Data Comp: Off Server Time: 2020-10-31 19:33:19 UTC
Valid HTML 5 and Valid CSS