[gmx-developers] Extracting coordinates from a trajectory
Florian Dommert
dommert at fias.uni-frankfurt.de
Thu Sep 4 10:59:56 CEST 2008
Hello,
since I try to create a new g_tool, with some help of Berk, many
thanks to him again, to calculate the dielectric constant for any kind
of systems, also ionic ones, some questions arise regarding the
coordinate transformation for unfolding the coordinates during the
analysis. My approach is more or less the same as in g_traj. If I want
to remove jumps, in other words unfold the trajectory, first the
remove_jump routine is used followed by rm_pbc, to obtain whole
molecules. After this procedure the COMs of the molecules are
calculated straightforward by multiplying the coordinates with the
corresponding factor and adding them up.
If only the averages of the translational (Md) and rotational (Mj)
part of the dipole moment are calculated from the folded coordinates,
in other words from the pure trajectory data, the procedure is
following. At first rm_pbc is used to make molecules whole, followed
by a straightforward calculation of the COMs of the molecules for the
translational part. At this stage pbc_dx provides the distance from
the COMs to the center of the box. This distance is multiplied by the
charge of the molecule and the results are added up.
Finally the sum over the positions of the COMs of the molecules times
their charge represents the translational fraction of the total dipole
moment of a simulation cell.
Now some problems arise. The tool also calculates the Mean square
displacement of Mj, that allows to determine the conductivity and <
Mj^2 >. However if I use unfolded coordinates to calculate the MSD and
compare the value of <Mj^2> extracted from the the MSD of Mj and from
the simulation snapshots, the difference is very large, regardless if
folded or unfolded coordinates are taken into respect for the
calculation. Although the simulated values almost correspond perfectly
to a straight line as expected, resulting in a very small error.
The next strange thing shows up, if folded coordinates are used to
calculate the MSD, what is physically senseless. In this case a curve
like a MSD is obtained, however it seems that the statistics are very
bad. A fit of a straight line gives me a < Mj^2 > that corresponds
almost with the value from the averaged data from the folded set, that
was also used for the MSD here.
Finally the trouble is, that <Mj^2> from the fit of the unfolded
coordinates yields an expected value for the dielectric constant,
whereas the simple average gives numbers two orders of magnitude larger.
Finally a questions arise. How far does the pressure coupling affect
the unfolding prodecure, if the size of the box is always changing ?
Can this problem be circumvented with the new mdp options refscale.
Unfortunately I was not able to find a manual describing the new
options. The cvs manual seems to be buggy. For example the
transformation of the dihedral potentials to RB-potentials was not
correct that last time I looked into this CVS module (about 4 month
ago).
Here is the code of g_current, that is already in the cvs. However
here is an improved version, that is able to do an Einstein-Helfand
fit and an calculated the conductivity from the current
autocorrelation function by numerical integration first, that is
enhanced by n analytic integral of a fitted function a*t^(-b). These
two option are not present in the cvs-Version. To this is ask somebody
with write access to put the code into the cvs.
At the end a last question regarding g_velacc. David posted why not to
include the calculation of the VACF of molecules. As I wanted to
implement this by just adding a routine preparing the corresponding
factor for the velocities, I realized that you get a segfault, if you
want to calculate any kind of correlation function for molecules,
since the atom indices are used instead of molecule indices. I do not
wont to create buggy code so my question. Is the routine atom2molindex
the correct way to solve the problem ?
Hopefully my post describes the problem in the right way, summarized I
am asking myself why do I get so large values from averaging ?
Best Regards,
Flo
-------------- next part --------------
A non-text attachment was scrubbed...
Name: gmx_current.c
Type: application/octet-stream
Size: 20769 bytes
Desc: not available
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-developers/attachments/20080904/1da8b90b/attachment.obj>
-------------- next part --------------
--
Florian Dommert
Dipl.-Phys.
Computational and Theoretical Softmatter & Biophysics group
Frankfurt Institute for Advanced Studies
Johann-Wolfgang-Goethe University
Ruth-Moufang-Str. 1
60438 Frankfurt am Main
Phone: +49(0)69 / 798 - 47522
Fax: +49(0)69 / 798 - 47611
EMail: dommert at fias.uni-frankfurt.de
Home: http://fias.uni-frankfurt.de/~simbio/Florian_Dommert
-------------- next part --------------
A non-text attachment was scrubbed...
Name: PGP.sig
Type: application/pgp-signature
Size: 194 bytes
Desc: This is a digitally signed message part
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-developers/attachments/20080904/1da8b90b/attachment.sig>
More information about the gromacs.org_gmx-developers
mailing list