I am a novice in MD simulation. I want to calculate solvation time correlation function Cs(t) by MD simulation but I don't know the procedure.Can anybody help me ?
What system are you interested in? Chromophore + solvent, protein, ...? The term 'solvation time correlation function (TCF) ' is usually used for a solute that undergoes some kind of excitation (change of partial charges, geometry or size), so the rest of my answer is concerned solely with according systems.
The basic procedure is: Set up the desired system and let it equilibrate at the desired temperature (this is very important, as changes due to nonequilibrated systems will spoil your results). Then, in your production run, save the coordinates at the desired time-intervalls. All that follows is post-processing of your trajectory, use for example python+MDAnalysis:
What is the solvation energy you are interested in? In polar systems, one usually uses the Coulomb energy. In nonpolar systems, one uses the Lennard-Jones energy. For every timestep, calculate the change in interaction energy Delta U (Coulomb, Lennard-Jones, or both) between solute and solvent in the ground and excited state (make sure you correctly treat period boundary conditions). As the time-dependence of the interaction energy arises solely from the change in coordinates of solute and solvent, this is easy to calculate if you have saved the coordinates of your trajectory (you also need the partial charges, and LJ-parameters). Calculate the mean of the change in interaction energy Delta U and subtract it. Now, calculate the time correlation function (correlate (Delta U - mean)(t) with (Delta U - mean (0)).
You should also consider whether TCFs describe a solvation response well, (this is, if linear response theory applies), otherwise your results will not be meaningful.
Thanks a lot for your comment.It will surely help me.But I have few confusions which I want to clear. First of all my system is chromophore+ polar solvent and I want to calculate LJ and Coulomb energy separately for my system.Is it possible ,then how? You are talking about mean of the change of interaction energy delta U, how can I calculate it?
Does here Delta U indicate total energy or potential energy (using ener.edr file which one I select total energy or potential energy)? And please also advise me how I do excited state simulation (all I have done it is of ground state ).
You are only interested in the interaction energy between chromophore and solvent, this is usually not given by the program itself, as usually only the total or potential energy of the whole system is printed. We program this by ourselves. E.g. for the Coulomb energy:
Delta U = sum_solvent_atoms sum_chromophore_atoms ((q_solvent_atom delta q_chomophore_atom)/r)
calculated for each time step separately
You end up with a large array of Delta U at all timesteps. For the mean of Delta U, simply sum up all array entries and divide by the number of entries.
Does your chomophore change in size? Otherwise the change in Lennard-Jones energy is zero and it is useless the calculate it.
If you are interested in time-correlation functions, you do not need to run excited state simulations. All you need to know is how the partial charges (and Lennard Jones parameter if you want to do that) change upon excitation.
We published a number of articles on that topic, you might find some help there as well (e.g. J Phys Chem (2016) 145, 164506, J Phys Chem (2016) 145, 164507, Phys Chem Chem Phys (2017) 19, 10940, Phys Chem Chem Phys (2018) 20, 5246). Also have a look on Mark Maroncellis early work (1990--), this will help you, too!