Amber-MD / cpptraj

Biomolecular simulation trajectory/data analysis.
Other
138 stars 64 forks source link

Export dataset based on maths and distances #1104

Closed H-EKE closed 1 month ago

H-EKE commented 1 month ago

Hi,

I would like to calculate with cpptraj an activation index which is based on a set of distances and coeffiecients

My idea was the follwoing script

parm prod1/structure.psf
for i=1;i<2;i++
trajin prod$i/output.xtc 1 250 1
strip :SOD,CLA,TIP3,TIP,POPC
autoimage
rms first F$i @CA 
dataset F$i legend Rep$i
# A100 index
# A100 = -14.43x(V1.53-L7.55)-7.62(D2.50-T.3.37)
#       + 9.11x(N3.42-I4.42)-6.32x(W5.66-A6.34)-5.22(L6.65-Y7.35)+278.88
# V1.53-L7.55
distance A$i :1088@CA :1370@CA
#D2.50-T.3.37
distance B$i :1116@CA :1157@CA
#N3.42-I4.42
distance C$i :1162@CA :1189@CA 
#W5.66-A6.34
distance D$i :1268@CA :1306@CA
#L6.65-Y7.35
distance E$i :1330@CA :1350@CA
create test.dat A100
A100 = -14.43*(A$i)-7.62*(B$i)+9.11*(C$i)-6.32*(D$i)-5.22*(E$i)+278.88 
go
clear trajin
done

but my dataset A100 is alway empty. Is there anyway to export it as dataset, so at the end I have somehting like this:

Frame | Replicate 1 | Replicate 2| Replicate 3 A100 value| A100 value | A100 value

Thanks in advance!

drroe commented 1 month ago

Hi,

So there are 2 issues here. The first is that the actions are actually queued up for the next trajectory processing run, but the calculation statement (i.e. A100 = ...) is executed immediately, so at the time of the calculation (before go the distance commands haven't yet executed and the corresponding arrays are empty. The second issue is that since you're doing this within a loop, you need to give the A100 result a different name each time (assuming you will increase the number of iterations in the future) or it will be overwritten.

I haven't tested this input, but something like this should work:

parm prod1/structure.psf
for i=1;i<2;i++
  trajin prod$i/output.xtc 1 250 1
  strip :SOD,CLA,TIP3,TIP,POPC
  autoimage
  rms first F$i @CA 
  dataset F$i legend Rep$i
  distance A$i :1088@CA :1370@CA
  distance B$i :1116@CA :1157@CA
  distance C$i :1162@CA :1189@CA 
  distance D$i :1268@CA :1306@CA
  distance E$i :1330@CA :1350@CA
  # Execute trajectory processing to actually calculate the distances
  go
  # Calculate the parameter from the distances
  A100.$i = -14.43*(A$i)-7.62*(B$i)+9.11*(C$i)-6.32*(D$i)-5.22*(E$i)+278.88
  clear trajin
done
writedata test.dat A100.*

Hope this helps!

PS - I know it can be a bit confusing how Actions and Analyses are queued for trajectory processing and other commands are executed immediately. This is legacy behavior from PTRAJ and I chose to follow it so as to remain as backwards compatible with PTRAJ as possible. Hopefully my response clears things up a bit.

H-EKE commented 1 month ago

Dear @drroe

Thanks a lot, it worked perfectly. You saved me once again!

Go and run are always quite confusing for me, speciallly when one has a loop before.