This code analyzes molecular dynamics trajectory files. It currently can read .xyz or .xtc files and works with water, methanol, and acetonitrile. It can easily be extended to accommodate any simulation of a pure polar liquid. This code was used to produce the results in this publication:
Elton, D. C. and Fernandez-Serra, M.-V. “The hydrogen bond network of water supports propagating optical phonon-like modes” Nat. Comm. 7, 10913 (2016)
The epskw.f90 code can compute the following quantities:
- static longitudinal and transverse nonlocal susceptibiliies, chi_L(k,0), chi_T (k,0)
- longitudinal and transverse polarization correlation functions phiL(k, t), phiT(k, t)
- distance decompositions chi_L(k,0,R), chi_T (k,0,R), phiL(k, t, R), phiT(k, t, R), for the smallest k-value accessible in the system.
- static structure factors S_L(k,0) and S_T(k,0) and (experimental feature, may not be normalized correctly)
The /kw/kw.m MATLAB code can compute:
- nonlocal susceptibiliies, chi_L(k,w), chi_T(k,w), using the output files from epskw.f90
- plots of tau vs k for exponential relaxation or tau & omega vs k for oscillator peaks
To compile you need the xtc library libxdrf.a. A precompiled version is included in this repo, but chances are it will not work on your system and you will need to compile your own version. To obtain this library, download the source code from ftp://ftp.gromacs.org/pub/contrib/xdrfile-1.1.tar.gz, edit the makefile as necessary and make. Copy libxdrf.a into the /epskw folder. Run "make" to compile
Setting the input options
The included file epskw.inp gives an example of an input file, which is fixed format. To use a small set of k values (10-20), set use small k set ? to .t.. To use a dense set of k-points, select a maximum number of points to use using the diagonals and edges of the box.
The following forcefield models are supported:
- TIP4P, TIP4P2005, TIP4P2005f, SPCE, TIP3P,
- TTM3F (requires two file containing charges of all atoms in each timestep in the order OHHOHHOHH.., and a file containing the electronic polarization vectors. no spaces between timesteps. Use the same prefix as the coordinate file, but drop the suffix ".xyz") and add "Edip.dat" and "chgs.dat" instead.)
- methanol (GAFF)
- methanolH1 (H1+3 parameters)
- acetonitrile (GAFF)
To add your own forcefield model, edit main_stuff.f90.
Running the program
to run, use:
./epskw.x < epskw.inp
Plotting the results
Output files are outputted in a form readible by xmgrace. To view the figures use:
xmgrace -nxy *output_filename*
Copyright 2014-2016 Daniel C. Elton
This software is licensed under The MIT License (MIT) Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.