[gmx-developers] Extracting coordinates from a trajectory

[Florian Dommert]

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>