Frequency calculation for QM/MM

[jokurian]

Hi admins,
I am trying to do some QM/MM dynamics using ORCA-Tinker in SHARC. For a normal dynamics simulation in SHARC, I can easily do a frequency calculation in any QM software and give the molden format to SHARC but for a QM/MM simulation, this step is difficult.
Do you have any suggestion on how to do this frequency calculation or does SHARC provide any code to do it?

Thank you,
Jo Sony Kurian

[maisebastian]

Dear Jo Sony Kurian,
what is your goal of doing a frequency calculation in QM/MM? Are you interested in vibrational spectroscopy? Or do you simply want to sample initial conditions?

In general, a "normal" frequency calculation is essentially based on the assumption that a single harmonic oscillator can describe all important parts of the potential energy surface (PES). This is typically true for small, rigid molecules without torsions, flexible side chains, etc. In this case, one can approximate the PES with a second-order Taylor series around a single minimum and analyze this simpler potential instead of the true PES.
However, in QM/MM systems, this is usually not a good idea. The presence of many independent molecules leads to a huge conformational freedom, giving many local minima in the PES. This means that you cannot choose a single minimum and expand your Taylor series around it. Furthermore, normal modes are linear modes, but the solvent diffusional modes cannot be described by linear coordinates. Hence, a normal mode frequency analysis (in the way normally done in quantum chemistry) typically does not make much sense for QM/MM systems.

Instead, in such flexible systems, everything is about sampling. Essentially, you run a long molecular dynamics trajectory of your QM/MM system, and that provides you with the relevant information about the PES. You can get the vibrational spectra by the Fourier transformation of the autocorrelation function of the dipole moment of your molecule. And you can get initial conditions by randomly picking snapshots from the MD trajectory.

Best regards,
Sebastian

[jokurian]

Hi Sebastian,
Thank you for the reply and the detailed description.
My aim is to sample initial conditions and doing md simulation seems to be the easy way. I see that SHARC provides "sharctraj_to_initconds.py" to sample from "SHARC trajectories" but I dont understand how to do a simple md calculation with SHARC without doing the entire TSH. I would appreciate it if you could give some idea about it.
Thank you,
Jo

[maisebastian]

Hi Jo,
there are actually two scripts intended to get initial conditions from MD simulations. On one side, you can use sharctraj_to_initconds.py, but this does not work in your case, as you need to have already finished a SHARC simulation before you can use the script. This script is intended to restart trajectories with other settings. For example, you could run trajectories in SHARC in only the ground state, then after 200fs collect the last snapshot from all trajectories and restart them in the S3 state.

For your application, you can have a look at amber_to_initconds.py. This allows you to setup an AMBER MD simulation, run a trajectory and then extract snapshots from this trajectory to use as initial conditions for SHARC. Here, the use of classical MD is strongly recommended, because only over many ps or ns is it possible to properly sample the solvent phase space. So, you could do:

1. setup your system for classical MD in AMBER. You will need force fields for all molecules, which might require that you parametrize something.
2. Construct your system (molecule plus solvent shell plus anything else)
3. Minimize, heat, equilibrate, run an MD trajectory
4. Use amber_to_initconds.py to convert the trajectory into an initconds file
5. Use the initconds file to setup QM/MM SHARC trajectories

This seems complicated, but for QM/MM there is really no other useful way of generating initial conditions than running a ground state MD simulation. And for this, only classical MD is fast enough to sample properly.

If you do not have AMBER, you can use your favorite other MD code. Unfortunately, we do not have converter scripts yet, but adapting amber_to_initconds.py should be doable with moderate effort. You need to have access to all coordinates and velocities at sufficiently spaced time intervals. Then it is just a matter of reformating an ASCII file.

Best,
Sebastian

[jokurian]

Dear Sebastian,
Thank you for the reply.
I am trying to convert the md data from an MD code to initconds. Is the energies given under the geometry and velocity necessary? Does the following calculations depend on any of those values?

Thank you,
Jo

[maisebastian]

Dear Jo,
Having the energy lines there is currently necessary so that the file can be correctly parsed, but the numerical values are only for information. You can safely set the numbers to zero. Once you read in the initconds file with excite.py, the Ekin will be updated automatically.

Best,
Sebastian

[jokurian]

Dear Sebastian,
Thank you for the help!

Closing the thread