[gmx-users] request for future tools

[Mark Abraham]

David Mobley wrote:
> Dear all, (especially developers),
> 
> I just wanted to put in a request for future development of g_energy and 
> similar tools:
> 
> For scripting purposes, it would be really useful if the output 
> categories were not dependent on the details of the run. For example, if 
> I run g_energy on one of my free energy calculations with orientational 
> and distance restraints, the potential enery is #15 and the dV/dlambda 
> is #20. But without my orientational and distance restraints, the 
> potential energy is #10 and the dV/dlambda is #15. It would be much 
> better if both cases could be consistent and, in the case without 
> orientational restraints, the values corresponding to restraint energies 
> would simply return zero.

I agree that it would be a good feature. It is "only" a matter of 
establishing a mapping of all possible output terms to an integer. This 
isn't easy if you have multiple (nonstandard) groups, because the 
interaction energies of the groups is reported in the .edr file and the 
number and identity of the groups isn't known at compile time. That can 
be gotten around by having their indices all greater than the largest 
index of a term known at compile time, or all letters, etc.

> At present, it's really obnoxious from the perspective of writing 
> analysis scripts, as my analysis scripts either need to know whether 
> distance restraints are turned on or not (in order to extract the 
> correct energies, if I want, say, the potential energy) or they need to 
> somehow run g_energy, look for the potential, and figure out what value 
> to input to get it. Does anyone have an easier option?

I don't know whether you would agree that the solution I suggest below 
is "easier", but it certainly works. I use the tool expect (found on all 
modern *NIXES) to do more advanced versions of the oft-recommended

echo "3 0" | g_energy -f output.edr

schemes, so that I can deal with a limited number of known .edr formats. 
Every time the script breaks you can add in a new format following the 
example below and use to your heart's content. The example below is an 
old attempt at extracting the total Coulomb energy from a gromacs 3.2.1 
.edr file that could have been calculated with various electrostatic 
methods. You could modify this suitably to treat your known range 
variation (email me if you want the file as an attachment). In theory an 
expect script could digest the g_energy output and synthesize the 
responses, but that would be overkill.

Regards,

Mark

-----------------

#!/usr/bin/expect -f

# Usage: expect this_script.tcl command_to_automate [args]

set force_conservative 0  ;# set to 1 to force conservative mode even if
			  ;# script wasn't run conservatively originally
if {$force_conservative} {
	set send_slow {1 .1}
	proc send {ignore arg} {
		sleep .1
		exp_send -s -- $arg
	}
}

set timeout -1

# now construct and run the sub-process
set cmd [list spawn ]
foreach a $argv {
   lappend cmd $a
}
eval $cmd

match_max 100000
while 1 {
     expect {
	"End your selection with 0\r
    1=         Angle   2=   Proper Dih.   3=Ryckaert-Bell.   4= 
LJ-14\r
    5=    Coulomb-14   6=       LJ (SR)   7=  Coulomb (SR)   8=  Coul. 
recip.\r
    9=     Potential  10=   Kinetic En.  11=  Total Energy  12= 
Temperature\r
   13=Pressure (bar)  14=        Vir-XX  15=        Vir-XY  16= 
Vir-XZ\r
   17=        Vir-YX  18=        Vir-YY  19=        Vir-YZ  20= 
Vir-ZX\r
   21=        Vir-ZY  22=        Vir-ZZ  23= Pres-XX (bar)  24= Pres-XY 
(bar)\r
   25= Pres-XZ (bar)  26= Pres-YX (bar)  27= Pres-YY (bar)  28= Pres-YZ 
(bar)\r
   29= Pres-ZX (bar)  30= Pres-ZY (bar)  31= Pres-ZZ (bar)  32= 
#Surf*SurfTen\r
   33=          Mu-X  34=          Mu-Y  35=          Mu-Z  36= 
T-Protein\r
   37=  Lamb-Protein\r" { send -- "5 7 8 0\r" }
	"End your selection with 0\r
    1=         Angle   2=   Proper Dih.   3=Ryckaert-Bell.   4= 
LJ-14\r
    5=    Coulomb-14   6=       LJ (SR)   7=  Coulomb (SR)   8= 
Potential\r
    9=   Kinetic En.  10=  Total Energy  11=   Temperature  12=Pressure 
(bar)\r
   13=        Vir-XX  14=        Vir-XY  15=        Vir-XZ  16= 
Vir-YX\r
   17=        Vir-YY  18=        Vir-YZ  19=        Vir-ZX  20= 
Vir-ZY\r
   21=        Vir-ZZ  22= Pres-XX (bar)  23= Pres-XY (bar)  24= Pres-XZ 
(bar)\r
   25= Pres-YX (bar)  26= Pres-YY (bar)  27= Pres-YZ (bar)  28= Pres-ZX 
(bar)\r
   29= Pres-ZY (bar)  30= Pres-ZZ (bar)  31= #Surf*SurfTen  32= 
  Mu-X\r
   33=          Mu-Y  34=          Mu-Z  35=     T-Protein  36= 
Lamb-Protein\r" { send -- "5 7 0\r" }
	"End your selection with 0\r
    1=       LJ (SR)   2=  Coulomb (SR)   3=     Potential   4= 
Kinetic En.\r
    5=  Total Energy   6=   Temperature   7=Pressure (bar)   8= 
Vir-XX\r
    9=        Vir-XY  10=        Vir-XZ  11=        Vir-YX  12= 
Vir-YY\r
   13=        Vir-YZ  14=        Vir-ZX  15=        Vir-ZY  16= 
Vir-ZZ\r
   17= Pres-XX (bar)  18= Pres-XY (bar)  19= Pres-XZ (bar)  20= Pres-YX 
(bar)\r
   21= Pres-YY (bar)  22= Pres-YZ (bar)  23= Pres-ZX (bar)  24= Pres-ZY 
(bar)\r
   25= Pres-ZZ (bar)  26= #Surf*SurfTen  27=          Mu-X  28= 
  Mu-Y\r
   29=          Mu-Z  30=     T-Protein  31=  Lamb-Protein\r" { send -- 
"2 0\r" }
	"End your selection with 0\r
    1=       LJ (SR)   2=  Coulomb (SR)   3=  Coul. recip.   4= 
Potential\r
    5=   Kinetic En.   6=  Total Energy   7=   Temperature   8=Pressure 
(bar)\r
    9=        Vir-XX  10=        Vir-XY  11=        Vir-XZ  12= 
Vir-YX\r
   13=        Vir-YY  14=        Vir-YZ  15=        Vir-ZX  16= 
Vir-ZY\r
   17=        Vir-ZZ  18= Pres-XX (bar)  19= Pres-XY (bar)  20= Pres-XZ 
(bar)\r
   21= Pres-YX (bar)  22= Pres-YY (bar)  23= Pres-YZ (bar)  24= Pres-ZX 
(bar)\r
   25= Pres-ZY (bar)  26= Pres-ZZ (bar)  27= #Surf*SurfTen  28= 
  Mu-X\r
   29=          Mu-Y  30=          Mu-Z  31=     T-Protein  32= 
Lamb-Protein\r" { send -- "2 3 0\r" }
	eof exit
     }
}