[gmx-users] Re: deshuf.ndx doesn't work for forces?

Chris Neale chris.neale at utoronto.ca
Fri Jan 25 00:48:08 CET 2008


>I think you need to add -force option.
>Like trjconv -f a.trr -force -o b.trr

Oh, yes. Thank you. I didn't know about that. Here is a script that I used to deshuffle a .trr file (It was never shuffled, this is just a
proof of concept). Indeed the trjconv that receives a .ndx file that is in a non-ascending order does change the x and v but not the f.

Using this reverse reordering it is simple to see if the correct modification has taken place.

Below I post 3 thins:

1) A diff of my source code modification to gmx_trjconv.c 3.3.2
2) Some parsed output from my final diff test
3) The code for a test

Please note that I do not in any way endorse the code change here. I did my best, but I don't work too much with the gromacs 
code so I could easily have broken something. I tried only to copy what is done with velocities.


####################################################
1) A diff of my source code modification to gmx_trjconv.c 3.3.2

$ diff gmx_trjconv_orig.c gmx_trjconv.c
663c663
<   rvec         *xmem=NULL,*vmem=NULL;
---
>   rvec         *xmem=NULL,*vmem=NULL,*fmem=NULL;
999a1000
>         snew(fmem,nout);
1205a1207
>               frout.f = fmem;
1208a1211
>                 if (bForce && fr.bF) copy_rvec(fr.f[index[i]],frout.f[i]);
 

####################################################
2) Some parsed output from my final diff test

>       f[    0]={-1.08121e+02, -1.03828e+03, -4.96866e+02}
>       f[    1]={-8.78225e+01, -3.77849e+02, -3.16504e+02}
>       f[    2]={ 6.70155e+01,  1.78264e+03,  7.23227e+02}
>       f[    3]={-7.33609e+02,  1.13756e+02,  2.29893e+02}
>       f[    4]={-4.07243e+02,  4.48937e+02,  1.64252e+02}
...
>       f[ 1913]={ 1.23753e+03,  2.47442e+02,  5.33197e+01}
>       f[ 1914]={-1.16692e+02, -5.33725e+01,  1.52119e+02}
>       f[ 1915]={-1.26504e+02, -3.05828e+02,  8.84993e+01}
>       f[ 1916]={ 4.86626e+01,  1.89065e+02,  1.42950e+02}
>       f[ 1917]={ 2.51086e+01, -1.58804e+02,  4.47615e+01}


<       f[    0]={ 2.51086e+01, -1.58804e+02,  4.47615e+01}
<       f[    1]={ 4.86626e+01,  1.89065e+02,  1.42950e+02}
<       f[    2]={-1.26504e+02, -3.05828e+02,  8.84993e+01}
<       f[    3]={-1.16692e+02, -5.33725e+01,  1.52119e+02}
<       f[    4]={ 1.23753e+03,  2.47442e+02,  5.33197e+01}
...
<       f[ 1913]={-4.07243e+02,  4.48937e+02,  1.64252e+02}
<       f[ 1914]={-7.33609e+02,  1.13756e+02,  2.29893e+02}
<       f[ 1915]={ 6.70155e+01,  1.78264e+03,  7.23227e+02}
<       f[ 1916]={-8.78225e+01, -3.77849e+02, -3.16504e+02}
<       f[ 1917]={-1.08121e+02, -1.03828e+03, -4.96866e+02}


####################################################

3) The code for a test


#!/bin/bash

grompp -f adp.mdp -c adp_md0.gro -p adp.top -n adp.ndx -o adp_md1.tpr
mdrun -deffnm adp_md1

#create a bogus .ndx file for no deshuffling
numAtoms=`head -2 adp_md0.gro | tail -1 | awk '{print $1}'`
{
echo '[ DeShuffle ]'
for((i=1,c=0; i<=numAtoms; i++)); do
  echo -n "$i "
  let "c=$c+1"
  if((c>10)); then
    c=0
    echo ""
  fi
done
}> deshuf.ndx

#create a bogus .ndx file for inverted deshuffling
{
echo '[ DeShuffle ]'
for((i=numAtoms,c=0; i>=1; i--)); do
  echo -n "$i "
  let "c=$c+1"
  if((c>10)); then
    c=0
    echo ""
  fi
done
}> deshuf2.ndx


trjconv -f adp_md1.trr -o adp_md1_deshuffled.trr -s adp_md1.tpr -n 
deshuf.ndx -force
gmxdump -f adp_md1.trr > dump.reg
gmxdump -f adp_md1_deshuffled.trr > dump.deshuffled
diff dump.reg dump.deshuffled > diff.regDeshuffled.dump

trjconv -f adp_md1.trr -o adp_md1_deshuffled2.trr -s adp_md1.tpr -n 
deshuf2.ndx -force
gmxdump -f adp_md1_deshuffled2.trr > dump.deshuffled2
diff dump.reg dump.deshuffled2 > diff.regDeshuffled2.dump




----------------------------------------
>/ Date: Thu, 24 Jan 2008 17:39:45 -0500
/>/ From: chris.neale at utoronto.ca <http://www.gromacs.org/mailman/listinfo/gmx-users>
/>/ To: gmx-users at gromacs.org <http://www.gromacs.org/mailman/listinfo/gmx-users>
/>/ Subject: [gmx-users] Re: deshuf.ndx doesn't work for forces?
/>/=20
/>/ My apologies, The commands should have been:
/>/=20
/>/ trjconv -f a.trr -o b.trr
/>/ gmxdump -f a.trr> a.dump
/>/ gmxdump -f b.trr> b.dump
/>/ diff a.dump b.dump> ab.diff
/>/=20
/>/ Chris Neale wrote:
/>>/ OK, I can confirm this for 3.3.1. However, I don't think that the=20
/>>/ problem is one of shuffling or sorting but simply that
/>>/ trjconv -f a.trr -o b.trr creates an exact copy of the .trr file=20
/>>/ *without* the forces. You can easily verify this by doing this:
/>>/
/>>/ trjconv -f a.trr -o b.trr
/>>/ gmxdump -f a.trr> a.dumo
/>>/ gmxdump -f b.trr> b.trr
/>>/ diff a.trr b.trr> ab.diff
/>>/
/>>/ Note that b.trr is of smaller size then b.trr and the ab.diff file=20
/>>/ clearly shows that it is the forces that are missing.
/>>/
/>>/ I actually have never tested that before now so it appears that my=20
/>>/ g_desort routine will lead to the loss of force information at each=20
/>>/ re-sorting event by nature of passing through trjconv -o .trr.
/>>/
/>>/ I did take a quick look at the gmx_trjconv code but I am unable to=20
/>>/ determine the particular place where the code might be changed.
/>>/ However, I think that I have isolated it to the need to add an output=20
/>>/ of force data from the -o .trr option to trjconv.
/>>/
/>>/ Chris
/>>/
/>>/
/>>/
/>>/ --- Original Message ---
/>>/
/>>/ Hi,
/>>/ I think after the "trjconv -f *.trr -n deshuf.ndx" operation,  the new=
/=20
>>/ or=3D
/>>/ ders of x and v in the trr file are correct.=3D20
/>>/ Here "correct" means the order is the same as the intitial *.gro file=20
/>>/ you=3D
/>>/ use for grompp command. If the output
/>>/ of trjconv is gro/pdb, you're right that it's only deshuffled, but not=
/=20
>>/ de=3D
/>>/ sorted. But since there's no force in the gro
/>>/ file, the bug i mentioned is only in the case when you choose trr as=20
/>>/ the =3D
/>>/ output of your deshuffle operation.
/>>/ Lanyuan Lu
/>>/ ----------------------------------------
/>>>/ / Date: Thu, 24 Jan 2008 12:28:47 -0500
/>>/ />/ From: chris.neale at utoronto.ca=20
/>>/=20
/>>/ />/ To: gmx-users at gromacs.org=20
/>>/=20
/>>/ />/ Subject: [gmx-users] deshuf.ndx doesn't work for forces?
/>>/ />/=3D20
/>>/ />>/ Message: 2
/>>/ />>/ Date: Wed, 23 Jan 2008 18:56:22 -0500
/>>/ />>/ From: LuLanyuan=3D20
/>>/ />>/ Subject: [gmx-users] deshuf.ndx doesn't work for forces?
/>>/ />>/ To:=3D20
/>>/ />>/ Message-ID:=3D20
/>>/ />>/ Content-Type: text/plain; charset=3D3D"gb2312"
/>>/ />>/
/>>/ />>/
/>>/ />>/ Hello All,
/>>/ />>/ I found sth like a bug regarding the deshuffle index file. I did=20
/>>/ a  =3D20
/>>/ />>/ simulation using multiple cpus and -sort -shuffle options. After =
/=20
>>/ =3D20
/>>/ />>/ that I tried to use the deshuf.ndx file to recover the original=20
/>>/ atom =3D20
/>>/ />>/  order for my trr file by trjconv -n deshuf.ndx ... command.
/>>/ />/=3D20
/>>/ />/ As I understand it, the deshuf.ndx file does not allow you to=20
/>>/ desort, =3D20
/>>/ />/ merely to deshuffle.
/>>/ />/=3D20
/>>/ />>/ But by checking the output trr using the gmxdump, I found the=20
/>>/ oders =3D20
/>>/ />>/ of x  and v were deshuffled while the f oder remained unchanged.=20
/>>/ In  =3D20
/>>/ />>/ another word the values of forces were not consistent with those=20
/>>/ of  =3D20
/>>/ />>/ positions and velocities in the new trr file.
/>>/ />>/ Is this a bug?
/>>/ />/=3D20
/>>/ />/ So the order of the positions and velocities has changed, but are=20
/>>/ they =3D
/>>/ /=3D20
/>>>/ / correct? I do not believe that gromacs 3.3.1 is capable of=20
/>>>/ desorting, =3D20
/>>/ />/ nor is it supported. I have uploaded a tool called g_desort to the=
/=20
>>/ =3D20
/>>/ />/ users contribution section. It works well and is tested, but it's =
/=3D20
>>/ />/ usage is complicated when you want to load in .trr files via gromp=
/p.
>>/ />/=3D20
/>>/ />>/ And is there anyway to get the correct oder for forces also?
/>>/ />/=3D20
/>>/ />/ I suggest that you try to reproduce this problem while using=20
/>>/ -shuffle =3D20
/>>/ />/ but not -sort.
/>>/ />/=3D20
/>>/ />>/ Thanks a lot.
/>>/ />>/ Lanyuan Lu/
/>>





More information about the gromacs.org_gmx-users mailing list