This is a guide for novices who try to apply the vibrational frequency maps to molecular dynamics simulations for the first time. Basically, vibrational frequency maps are used to obtain a trajectory of physical quantities used to calculate the linear and nonlinear vibrational spectra from a MD trajectory. 

[1] First, you should have completed a MD simulation session to obtain a trajectory file of your system. There are many available MD simulation packages for simulating chemical system of your interest such as AMBER, CHARMM, NAMD, LAMMPS, Gromacs, and so on. I should assume that you have previous experience of running MD simulations or have a basic knowledge of MD methodology. If you are not familiar with MD simulations, please consult relevant textbooks, for example, “Understanding molecular simulation” written by Daan Frenkel and Berend Smit. The following is an example of input file for AMBER MD simulation package to run a simulation to obtain the trajectory:

&cntrl

  irest = 0, ntx = 5,

  nstlim = 1000000, nsnb = 10, dt = 0.001

  tempi= 0.0, temp0 = 300.0,

  ntt = 3, gamma_ln = 1.0,

  ntb = 1, ntp = 0,

  ntf = 2, ntc = 2, cut = 10.0,

  tol = 1.E-6,

  ntwx= 10, iwrap = 1,

 &end

This input file implements a simulation under the condition of a constant temperature, constant volume, constant number of particles with temperature coupling. One notable thing for running MD simulation is determining by what interval the snapshot structures are written to the trajectory file. In the above example, ‘ntwx=2’ means that the snapshot structures are written to the trajectory file every tenth step which means that the interval between the snapshot structures is 10 fs because the time step is 1 fs in this simulation. This can be a short interval and can be increased to reduce the size of the trajectory file. But too large interval should be avoided when the vibration frequency is high, since it might not be sufficient to sample the vibrational motion of the molecules.

[2] When you got the trajectory file of your system, then you should calculate the trajectory of the vibrational frequencies and other relevant quantities from your MD trajectory. It is in this step that the vibrational frequency maps are necessary. Let us assume, for example, that you have a MD trajectory of a N-methylacetamide (NMA) molecule in water and you are going to calculate the vibration spectrum of the amide I mode of the NMA molecule. Firstly, go to the frequencymap.org site and select the ‘Download’ tab. Type ‘amide’ in the search box. Then, you got a list of vibrational frequency map files for amide I mode for various cases such as NMA, peptide, proteins. Under the titles of the listed maps, there are brief descriptions. We choose the ‘amide I’ map file with the description “Amide I mode vibration of NMA molecule” having the file named “nma_ham_jcp_2003_118_34915.vbm”. The name of the file contains information of the molecule (nma) and the name of the first author (ham) of the paper published in the Journal of Chemical Physics, volume 118, page 34914 in 2003, which contains the contents related to the map parameters in this vibrational frequency map file. If you click ‘Download’ button on the left of the name of this vbm file, the vibrational frequency map file is downloaded.

[3] The downloaded vbm file contains reference paper in the “%reference” section. Reading of this paper before proceeding is highly recommended. The detailed format of the vbm file is given in the ‘Tutorial’ tab of the frequencymap.org site. Essentially, the ‘%numbers’ section says that the molecule has 12 atoms and 6 interaction sites. The ‘%structure’ section gives the list of the constituent atoms with their coordinates. The coordinates are given to show how the atoms in the molecule are indexed. The ‘%sites on’ section says that the 1^st, 5^th, 6^th, 7^th, 8^th, and 9^th atoms in the ‘%structure’ section are the interaction sites. The parameters for the interaction sites are given in the ‘%map param’ section.

[4] For each snapshot of your MD trajectory, calculate the electrostatic potential on the interaction sites of the NMA molecule due to the water molecules in the simulation box and obtain the amide I mode frequency by multiplying the parameters in ‘%map param’ section to the calculated potential values on the interaction sites of the NMA molecule. An example code for calculating the vibration frequency shifts for each snapshot of the MD trajectory using the parameters from the vbm file in python language might look like the following:

 

for site in range(number_of_site):

    potential[site] = 0

    for iwater in range(number_of_water_molecules):

        for iatom in range(number_of_atoms_in_a_water_molecule):

            potential[site] += charge[iatom][iwater]/abs(r[site]-r[iatom][iwater])

frequency_shift = 0

for isite in range(number_of_site):

    frequency_shift += potential[site]*parameter[site]

 

In the above code you should include the minimum image convention if your simulation uses periodic simulation box.

[5] If you successfully completed the step [4], now you have a time series of the amide I mode vibration frequency of the NMA molecule. You can use this data to calculate the linear and nonlinear vibration spectrum for the amide I mode vibration of NMA. The problem of how to calculate the linear and nonlinear vibration spectrum is beyond this guide. The recommended way to calculate the vibration spectrum is to use already established program packages. I would like to suggest the NISE (numerical integration of Schrodinger equation) program developed by Jansen group. More information can be found in https://github.com/GHlacour/NISE_2017.