Do you want to reach potential energy minimimum (via energy miminization) or equilibrate the system at a specific temperature (via molecular dynamics (MD) simulation) for further analysis (production MD)?

If you want to do the latter, there is no need for energy minimization.

@Martin, thanks...Although I tried once to skip the minimization step and directly start the equilibration but failed because of some issues. I will look to it again and it is a good idea if possible, because i want the actual positions of Zn ions for the equillibration step. I just want to equilibrate the system in different pH conditions. Therefore I am altering the protonation step to mimic the system.

Start the equilibration MD simulation (directly without energy minimization) at a low temperature (e.g., 5 K) and restrain the solute (including Zn ions). Gradually increase the temperature and decrease (and remove) the restraints.

---

In general, modeling of metal ions can be tricky (big charge; metal binding).

Also, the position of Zn ions can be affected by conformational changes of the Zn-binding amino acid residues so I would check that as well, especially early in the equilibration MD.

In some cases, water takes part in the binding of Zn ions; if this is your case, make sure to include the Zn-coordinating ordered water molecule(s) in your system.

Why would one skip energy minimization? A suitable restraint can easily be applied to maintain ion coordination. Taking a crystal structure and dunking it into a highly ordered box of water (as is standard procedure in most software) is very artificial and leads to clashes. EM is a prerequisite for just about any simulation. I see no reason to skip it.

@Martin. 

(e.g., 5 K) and restrain the solute (including Zn ions)...

How to restrain, what changes should I introduce in the configuration file to restrain them? 

Gradually increase the temperature and decrease (and remove) the restraints....

You mean I should run it for some time then stop the process and re-submit it again with a bit higher temperature and so on or is there any step which can automatically increase the temperature with out stopping the program. I never did this step before of starting with a low temperature and increase gradually therefore.

In some cases, water takes part in the binding of Zn ions...

Yes in my case Zn coordinates to OH ion also, though i have removed it but will also try to keep it now. 

@JUSTIN

You suggested a suitable restrain and EM is a prerequisite for any simulation. Is there any way to perform EM in NAMD or you are just asking me to try this whole new bunch of things. Sorry If i am mistaken.

In your original question, you stated that you did energy minimization but the structure changed quite a bit. This is because force field parameters for species interacting with divalent cations like Zn2+ are notoriously bad because the additive force field approximation is lacking in this case (no polarization or charge-transfer effects, which are quite significant in a real system).

There are numerous papers with reparametrizations of amino acid side chains and the ions themselves to try to overcome this, but these solutions are often highly specific to the exact protein/active site for which they were created, and thus are generally (in my mind) undesirable. They "fix" a problem but do not provide a real solution.

A restrained minimization is of course possible, in NAMD via the extraBonds keyword and suitable file, in which you define distance restraints between the Zn2+ and their ligating atoms. This is quite straightforward and is routinely done. I do not understand why Martin is suggesting that you skip EM altogether and would be curious to hear his thoughts. I think the approach he suggests might work but could possibly be prone to artificial forces that arise from an un-minimized structure.

Justin: Of course, if there are obvious clashes that would blow up the system in MD simulation then it makes sense to do some energy minimization (EM) first. But if one can do without EM then there is no need for EM.

---

I have done many hundreds of equilibration MD simulations with various systems (protein/DNA; DNA); all without EM.

I want to cook some meat. I happen to have some meat in my refrigerator. I open the refrigerator, grab the meat and --- put it in the freezer?

Ashfaq: I do not use NAMD (I use Q) but I am sure it is possible to set up equilibration MD simulation in NAMD as consisting of several short consecutive MD simulations, in which one can define the reference temperature, the temperature control, the restraints, . . .

@ Justin and Martin i searched NAMD user guide and found that if I set a constrain for certain atoms then those will be treated as fixed in equilibration state also. Here, my aim is to check the distance between two Zn ions in my simulation. If my understandings from the user guide are correct then for sure I do not need setting of constrains for Zn ions. So i have only one option available to treat my protein directly in the equilibration, although there are chances that high energy can affect the simulation. 

Next, I am thinking to minimize the protein alone then alter the minimized PSF and PDB file manually and accommodate those Zinc ions. I am not sure will it work or not but let me try.  But if worked, in that case i hope will not face any internal heating issue. 

Thanks to you both, i learnt too many things in this discussion. 

Hi Ashfaq,

I read everyone's reply and I am not able to find out why you don't expect any variation in relative position during em. I mean even you told that you want to look for the distance between two zinc ions and that means it is going to change and you want to look its dynamics. For any force field, when you perform unrestrained simulation, you always observe a deviation which then converges after some time. The convergence time varies from system to system and depends on the constituent, size, precision etc. This will happen to your system also and that's why you call it dynamics.

Even if you minimize the protein separately and then try to put zn back into the pocket through structure alignment or whatever method, you will just find that the position of zn is now perfectly in alignment with the crystals structure's position but it is not as the relative orientation of all the amino acid side chain around it had been changed. Watch your steps when you move ahead because it could be more dangerous as minimization in presence and absence  of Zn ions will be completely different. When you introduce Zn back into the pocket, you may end with a system which is not minimized.

What important is, you can heat your system slowly during equilibration, and can try to remove the restrain applied during equiliberation slowly and in steps not all of a sudden. This will prevent the effect of velocity of solvent molecules on protein and zn. You can also take zn ions with ordered water around it but I would suggest not not to take if any direct interaction of zn ion with residues is reported in crystal structure.

Thanks

Dear Kumar, I have three same structures solved in different pH conditions i.e. 5.5, 6.5 and 7.5. All the protein structures are similar in every aspects but only the distances between 2 Zinc ions are changed. I wanted to know whether this pH has a role, if yes then why not to check all these three in mimic different pH conditions. For pH 5.5, i did one change and that is protonation of both Histidine residues, while for 6.5 and 7.5 i am following one protonation for His. 

Now what I observed during minimization, the two Zinc ions in every structures moved away from each other. I am a bit confused here, because i was expecting that minimization will not treat Zinc that way and i will be able to equilibrate all the proteins in a so called native states.

Some one suggested me that try constrain minimization, therefore I put my question here. So basically i want to know the behavior of Zn ions in simulation.

Histidines ligating Zn2+ ions are going to have very noncanonical pKa values. Have you done a proper pKa calculation in the presence of Zn2+, or are you assuming that the pKa will be 6.5? Divalent metal ions are Lewis acids and will cause ligating residues to have a strongly down-shifted pKa. If you're protonating the histidines, then there's no amount of restraining or fiddling that will give you an unperturbed binding site.

Hi Ashfaq. It is very important to take care of H+ concentration and protonation states and transfer of protons among side chain of amino acids. In CHARM force field as far as I know, there exist a method which takes care of transfer of protons which is very important in constant pH simulations.

Thanks

This is not a feature of the CHARMM force field, but the CHARMM software can carry out constant-pH simulations. NAMD does not, as far as I am aware. This is a very different case than what is being attempted here.

Justin, No i did not calculate the pKa value, i just guessed it from the pH we used and assume 6.5 pKa at physiological state. Second, it seems that NAMD is not going to help me as you said that it can not carry constant pH simulations. Now the situation became more tricky and I guess out of my reach :-(. I hope this discussion will help others for sure in some extent. Thank you all for your time. 

The discussion here has gotten pretty far off-topic from where you started. A constant-pH simulation is one in which the protonation states of the titratable residues can vary over time, and titration coordinates are integrated according to a specified pH. You don't need to do this, necessarily. NAMD is perfectly capable of carrying out a fixed-topology simulation as you want to. You just need to calculate the appropriate pKa values of your protein's residues. You cannot (and should not!) assume that everything has a canonical value.